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UNCERTAINTY  ANALYSIS  FOR  FLUID  MECHANICS  WITH  APPLICATIONS 


ROBERT  W.  WALTERS*  AND  LUC  HUYSEt 

Abstract.  This  paper  reviews  uncertainty  analysis  methods  and  their  application  to  fundamental 
problems  in  fluid  dynamics.  Probabilistic  (Monte-Carlo,  Moment  methods,  Polynomial  Chaos)  and  non- 
probabilistic  methods  (Interval  Analysis,  Propagation  of  error  using  sensitivity  derivatives)  are  described 
and  implemented.  Results  are  presented  for  a  model  convection  equation  with  a  source  term,  a  model 
non-linear  convection-diffusion  equation;  supersonic  flow  over  wedges,  expansion  corners,  and  an  airfoil;  and 
two-dimensional  laminar  boundary  layer  flow. 

Key  words,  stochastic,  probabilistic,  uncertainty,  error 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction  and  Motivation.  In  the  past  few  years,  there  has  been  numerous  papers  appearing 
in  the  Computational  Fluid  Dynamics  (CFD)  literature  addressing  the  subject  of  credible  CFD  simulations 
(see  [1],  [2],  [6],  [7],  [10],  [21],  [30],  [33],  [35],  [36],  [38],  [39],  [41],  [42],  [49]).  In  fact,  the  May  1998  issue  of 
the  AIAA  Journal  devoted  a  special  section  on  the  topic.  Among  the  key  issues  discussed  in  that  section 
were:  Code  Verification,  Code  Validation,  Certification  and  Sources  of  Uncertainty.  One  of  the  primary 
reasons  for  the  increased  interest  in  uncertainty  management  is  its  application  in  risk-based  design  methods. 
The  structures  community  and  dynamics  and  controls  discipline  have  a  long  history  in  uncertainty  analysis 
whereas  the  computational  fluid  dynamics  community  is  a  newcomer  in  this  area  due  in  part  to  the  relative 
infancy  of  the  discipline  and  in  part  to  the  large  cost  of  CFD  simulations.  However,  with  the  increase  in 
computing  power  and  software  improvements  over  the  past  two  decades,  stochastic  CFD  is  coming  of  age. 
Consequently,  the  primary  purpose  of  this  paper  is  to  serve  as  an  introductory  guide  for  engineers  with  an 
interest  in  fluid  dynamics  applications  of  uncertainty  analysis  methods. 

Before  proceeding  to  review  the  basic  methods,  some  discussion  on  nomenclature  is  warranted.  In  this 
paper  we  adopt  the  AIAA  definitions  for  error  and  uncertainty  [3] ,  namely: 

Definition  1.1.  Error;  A  recognizable  deficiency  in  any  phase  or  activity  of  modeling  and  simulation 
that  is  not  due  to  lack  of  knowledge. 

Definition  1.2.  Uncertainty:  A  potential  deficiency  in  any  phase  or  activity  of  the  modeling  process 
that  is  due  to  lack  of  knowledge. 

These  definitions  recognize  the  deterministic  nature  of  error  and  the  stochastic  nature  of  uncertainty. 
Uncertainty  can  be  further  categorized  into  aleatoric  (or  inherent  uncertainty)  and  epistemic  (or  model) 
uncertainty.  Further  categorizations  are  possible  and  discussed  in  [29].  This  report  focuses  on  methods 
for  describing  and  propagating  parameter  uncertainty  in  models.  In  a  companion  report  [29],  we  describe 
methods  for  dealing  with  model  form  and  boundary  condition  uncertainty  for  the  non-linear  Burgers  equation. 
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under  NASA  Contract  No.  NASl-97046  while  the  author  was  in  residence  at  ICASE,  NASA  Langley  Research  Center,  Hampton, 
VA  23681-2199. 
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Table  1.1 


Source  of  Uncertainty  and  Error  in  CFD  Simulations  -  summarized  from  Oberkampf  and  Blottner,  Ref.  [35] 


Source 

Examples 

Physical  Modeling 

(Assumptions  in  the  PDF) 

Inviscid  Flow 

Viscous  Flow 

Incompressible  Flow 

Chemically  Reacting  Gas 
Transitional/Turbulent  Flow 

Auxiliary  Physical  Models 

Equation  of  State 
Thermodynamic  properties 
Transport  properties 

Chemical  models,  reactions,  and  rates 

Turbulence  model 

Boundary  Conditions 

Wall,  e.g.  roughness 

Open,  e.g.  far-field 

Free  Surface 

Geometry  Representation 

Discretization  &  Solntion 

Truncation  error  -  spatial  and  temporal 
Iterative  convergence  -  steady  state 
Iterative  convergence  -  time  dependent 
Geometry  Representation 

Ronnd-Off  Error 

Finite  -  precision  arithmetic 

Programming  &  User  Error 

In  some  instances  in  this  report,  we  have  used  a  more  general  definition  of  uncertainty  that  includes  error 
where  it  is  clear  that  no  ambiguity  arises. 

In  the  AIAA  special  issue,  Oberkampf  and  Blottner  [36]  group  sources  of  uncertainty  and  error  arising 
from  the  simulation  of  physical  phenomena  governed  by  PDE’s  into  four  broad  categories: 

1.  Physical  modeling 

2.  Discretization  and  solution  errors 

3.  Computer  round-off  error 

4.  Programming  errors 

An  examination  of  Table  1.1  shows  that  many  sources  of  error  and  uncertainty  arise  in  CFD  simulations. 
It  is  generally  believed  that  discretization  error,  geometric  uncertainty  and  turbulence  model  uncertainty  are 
the  largest  sources  of  uncertainty  in  modern  Reynolds-Averaged  Navier-Stokes  simulations  and  collectively 
account  for  much  of  the  scatter  observed  between  experimental  and  computational  data  [2].  Discretization 
error  has  been  studied  extensively  and  a  number  of  techniques  have  been  proposed  for  modeling  this  error 
(see  e.g.  [40],  [41],  [43]).  Grid  adaptation  schemes  frequently  use  these  models  for  improving  the  base 
grid.  The  impact  of  geometric  uncertainty  has  been  studied  by  Darmofal  [13]  for  compressor  blades  using 
probabilistic  methods.  Although  relatively  little  has  been  done  for  quantifying  turbulence  model  uncertainty, 
there  has  been  some  work  in  this  area  [8],  [20].  More  recently,  Godfrey  [19]  used  the  continuous  sensitivity 
equation  approach  to  rank  the  relative  importance  of  the  closure  coefficients  in  three  turbulence  models: 
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the  Baldwin-Lomax  algebraic  model,  the  Spalart-Allmaras  one-equation  model  and  the  Wilcox  two-equation 
k  —  Lu  turbulence  model. 

With  the  present  state  of  computational  resources,  some  sources  of  error  can  be  made  negligible.  For 
example,  for  essentially  all  one-  and  two-dimensional  steady,  inviscid  or  laminar  flows,  a  user  typically  has 
sufficient  hardware  to  drive  the  discretization  error  to  very  low  levels,  essentially  zero,  leaving  model  un¬ 
certainty  as  the  only  significant  source  of  uncertainty.  Over  time,  this  trend  will  continue  and  eventually 
encompass  a  large  class  of  three-dimensional  simulations.  However,  without  further  research  effort  concen¬ 
trated  on  uncertainty  estimation  and  management,  model  uncertainty  will  likely  remain  relatively  constant 
over  time  since  it  is  not  a  known  error  that  can  simply  be  reduced  with  additional  computing  power. 

In  the  next  section,  we  review  methods  that  can  be  used  in  CFD  simulations  for  dealing  with  error  and 
uncertainty.  Next,  a  range  of  problems  starting  with  some  simple  model  problems  and  ending  with  laminar 
boundary  layer  flow  are  presented. 

2.  Review  of  Uncertainty  Analysis  Methods.  In  this  section,  we  briefly  review  deterministic  and 
probabilistic  methods  for  uncertainty  analysis.  Two  deterministic  uncertainty  analysis  methods  -  1)  Interval 
Analysis  and  2)  Propagation  of  error  using  sensitivity  derivatives  and  three  probabilistic  methods  1)  Monte 
Carlo,  2)  Moment  methods,  and  3)  Polynomial  Chaos  are  summarized  below. 

2.1.  Interval  Analysis.  The  basic  idea  in  interval  analysis  is  to  perform  operations  on  input  intervals 
that  contain  the  set  of  all  possible  values  of  the  input  in  such  a  way  that  the  output  interval  consists  of  all 
possible  values  of  the  result  of  the  operations  performed  on  the  input.  Consequently,  interval  results  represent 
maximal  error  bounds  (i.e.,  worst  case  results).  One  of  the  most  appealing  things  about  Interval  Analysis 
(lA)  is  that  it  can  be  implemented  in  a  systematic  way  on  modern  computing  systems  such  that  the  details 
of  the  interval  operations  are  transparent  to  the  user.  Thus,  one  can  take  an  existing  simulation  tool,  such  as 
a  CFD  code,  and  immediately  implement  interval  analysis  provided  that  it  is  supported  by  the  programming 
environment.  However,  it  should  be  pointed  out  that  different  expressions  for  computing  an  interval  output 
quantity  can  result  in  different  interval  widths  even  if  the  expressions  are  mathematically  equivalent  for 
pointwise  input.  To  demonstrate  this,  consider  the  following  two  expressions  that  are  equivalent  for  point 
values, 

(2.1)  f{x)  =  and  g{x)  = 

1+x  1+^ 

Table  2.1  shows  the  results  of  performing  interval  analysis  for  these  two  functions  for  input  intervals,  x, 
defined  in  terms  of  the  interval  midpoint  value,  x,  and  uncertainty,  £,  by 

(2.2)  X  =  x[l  —  s,\  +  s]. 

Values  in  the  table  correspond  to  £  =  1/10.  Note  that  the  interval  width  associated  with  the  second  relation, 
g{x)^  is  substantially  smaller  than  the  width  found  by  evaluation  of  f{x).  This  shows  that,  although  existing 
software  can  take  immediate  advantage  of  interval  analysis,  significant  improvement  in  the  size  of  the  error 
bounds  may  be  possible  by  careful  design  and  construction  of  the  operations  within  the  software.  However, 
this  investment  may  not  be  prudent  since  probabilistic  methods  provide  much  more  information  than  can 
be  obtained  through  interval  analysis.  This  example  also  illustrates  one  other  point:  the  fact  that  different 
interval  results  can  be  obtained  is  not  related  to  the  precision  of  the  calculation.  Here,  rational  numbers  were 
used  to  carry  out  the  operations  exactly.  Different  interval  output  occurs  because  set  theoretical  operations 
-  union,  intersection,  complement  are  affected  by  the  structural  form  of  the  operations. 
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Table  2.1 

Interval  analysis  results  for  two  simple  expressions  that  are  equivalent  for  point  values  but  not  for  intervals. 


X 

X 

f(x) 

g(x) 

f(x) 

g(x) 

l|f(x)|| 

llg(^)l 

D 

1 

1 

■  9  11' 

.10'  10. 

] 

1 

2 

1 

2 

■3  11' 

.7’  19. 

•9  11' 

.19’  21. 

3 

5 

4 

'9  11' 
.8’  8  . 

5 

9 

5 

9 

■  9  11' 

.19’  17. 

■  9  11' 

.17’  19. 

7 

2 

3 

2 

1 

'27  33' 
.20^  20. 

] 

3 

5 

3 

5 

'27  33' 
.53’  47. 

'27  33' 
.47’  53. 

4 

7 

4 

1 

'63  77' 
.  40  ^  40 . 

] 

7 

11 

7 

11 

[ 

7  77 

13’  103 

] 

[ 

63  77 

103’  117 

■] 

9 

2 

2 

'9  11' 
.5’  5  . 

2 

3 

2 

3 

■9  11' 

.16^  14. 

■9  11' 

.14’  16. 
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Moreover,  interval  arithmetic  inside  iteration  loops  results  in  error  growth  each  iteration  without  some 
modification  to  the  base  algorithm  (see  Section  3).  Since  many  fiuid  dynamics  codes  rely  on  iteration  for  the 
solution  process,  this  further  detracts  from  the  use  of  this  approach.  In  a  probabilistic  modeling,  varying 
degrees  of  correlation  between  random  variables  can  readily  be  taken  into  account  in  the  analysis.  The 
practical  examples  of  Section  5.2  and  6.2  will  illustrate  this  point  and  indicate  that,  depending  on  the 
value  of  the  correlation,  the  uncertainty  may  be  larger  or  smaller  than  in  the  case  where  all  variables  are 
independent.  Since  interval  analysis  is  a  deterministic  method,  it  cannot  take  advantage  of  this  information 
and  must  necessarily  compute  the  widest  bounds. 

2.2.  Propagation  of  Error  using  Sensitivity  Derivatives.  Error  propagation  using  sensitivity 
derivatives  has  been  in  use  for  many  years  (see  e.g.  Dahlquist  and  Bjorck  [9]).  li  u  =  where 

is  the  independent  variable  with  error  associated  with  it,  then  a  deterministic  approximation  to  the 
error,  Au,  is  given  by 


A  computational  fiuid  dynamics  example  of  this  technique  is  presented  in  the  work  of  Turgeon  et.  ah,  [11],  in 
which  the  laminar  flow  of  corn  syrup  was  analyzed  subject  to  uncertainty  in  the  viscosity  model,  the  thermal 
conductivity  model,  the  geometry,  and  the  boundary  conditions.  In  their  work,  they  used  the  continuous 
sensitivity  equation  method  (see  Section  4)  to  evaluate  the  sensitivity  derivatives,  ■^,  in  Eq.(2.3).  This 
results  in  an  error  estimate,  Au,  which  in  their  work  was  shown  to  bound  the  experimental  data.  We  wish 
to  emphasize  that  this  approach  is  based  on  the  assumption  that  a  particular  input  interval  (for  example, 
one  of  their  inputs  was  an  experimentally  measured  temperature,  T  ±  2%),  contains  the  entire  uncertainty 
interval  due  to  this  variable.  We  will  contrast  this  with  a  probabilistic  framework  in  Section  2.4. 

2.3.  Monte  Carlo.  Although  there  are  many  different  applications  of  Monte  Carlo  simulation,  both 
deterministic  and  probabilistic,  the  focus  in  this  report  is  on  probabilistic  simulation  methods.  A  briefly 
history  of  the  method  is  given  by  Hammersley  and  Handscomb  in  [24]  and  summarized  here.  Arguably, 
the  development  of  the  method  and  name,  Monte  Carlo,  is  considered  to  have  its  start  around  1944  when 
von  Neumann  and  Ulam  performed  direct  simulation  of  neutron  transport  problems,  primarily  as  a  tool  for 
atomic  bomb  research.  Shortly  after  the  end  of  World  War  II,  Fermi,  Ulam,  von  Neumann  and  others  began 
applying  Monte  Carlo  methods  to  deterministic  problems.  By  1948,  Monte  Carlo  estimates  of  the  eigenvalues 
of  the  Schrodinger  equation  had  been  obtained.  Sometime  thereafter,  this  work  came  to  the  attention  of  Dr. 
Stephen  Brush  at  the  Livermore  Laboratory  who  subsequently  unearthed  a  1901  paper  by  Lord  Kelvin  [31] 
in  which  remarkably  modern  Monte  Carlo  techniques  appear  in  connection  with  the  Boltzmann  equation. 
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Fig.  2.1.  Typical  histograms  obtained  by  sampling  from  a  Normal  distribtution  with  a  mean  of  0.25  and  a  standard 
deviation  of  0.025  corresponding  to  a  coefficient  of  variation,  CoV  =  10%.  Left  -  1,000  samples.  Bight  10,000  samples. 

Apparently,  the  methods  were  obvious  to  Lord  Kelvin  and  consequently  his  focus  was  on  the  results.  Even 
prior  to  this  application,  there  are  isolated  accounts  of  the  method  [23] . 

Here  we  demonstrate  one  of  the  simplest  of  all  Monte  Carlo  methods,  referred  to  as  crude  (or  basic) 
Monte  Carlo.  In  this  approach,  the  basic  procedure  is: 

1.  sample  input  random  variable(s)  from  their  known  or  assumed  (joint-)  probability  density  function 

2.  compute  deterministic  output  for  each  sampled  input  value(s) 

3.  determine  statistics  of  the  output  distribution,  e.g.,  mean,  variance,  skewness,  ... 

The  statistics  of  a  distribution  (mean,  variance,  skewness,  kurtosis,  ...)  can  be  determined  from  the 
definition  of  the  expected  value  of  a  function  of  a  random  variable,  say  g{^),  namely 

(2.4)  JE[giO]  =  J  giOviOde 

where  p(^)  is  the  probability  density  function  of  the  distribution  that  describes  some  event  or  process  and 
the  integration  is  over  the  support  of  the  PDF.  The  mean  of  the  probability  distribution  (also  referred  to  as 
the  first  moment  about  the  origin)  is 

(2.5)  f  =  ]E[£]  =  J  ipiOde 

The  moment  about  the  mean  is  given  by 

(2.6)  lE[(^-f)1  =  Ji^-lYpiOde 

The  variance,  skewness,  and  kurtosis  are  related  to  the  2”'^,  3’''^,  and  4*^  moments  about  the  mean.  In 
some  cases,  the  integrals  can  be  evaluated  analytically,  in  others,  the  integrals  are  replaced  by  discrete  sums. 
In  the  application  problems  to  follow,  we  sampled  from  a  Normal  (Gaussian)  distribution  with  a  mean  p 
and  standard  deviation  rr.  The  probability  density  function  (PDF)  of  the  Normal  distribution,  pn{0j 
denoted  N[p,  a],  is  given  by: 


(2.7)  pn{0  =  ^  ■ 

VZttit 

Two  typical  samples  are  shown  in  Figure  2.1  in  which  the  Gaussian  shape  of  the  underlying  PDF  is  evident. 
Frequently,  it  is  convenient  to  use  the  standard  normal  variable,  A^[0, 1],  i.e.,  a  Gaussian  with  a  mean  of  0 
and  a  variance  of  1,  its  definition  follows  directly  from  Eq.  (2.7). 
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The  Monte  Carlo  method  has  the  property  that  it  converges  to  the  exact  stochastic  solution  as  the 
number  of  samples,  n  — >  oo.  However,  convergence  of  the  mean  error  estimate  is  relatively  slow  since  the 
standard  deviation  of  the  mean  scales  inversely  with  the  square  root  of  the  the  number  of  samples. 


(2.8) 


\Ai 


Substantial  efficiency  improvements  over  the  basic  scheme,  known  as  variance  reduction  techniques,  are 
reported  in  the  literature  (see  e.g.  [22],  [24],  [25],  [32]).  The  two  primary  mechanisms  for  improving  the 
basic  procedure  are  importance  sampling  and  correlation  methods. 


2.4.  Moment  Methods.  A  number  of  applications  using  moment  methods  have  appeared  in  the  litera¬ 
ture  involving  CFD  simulations  (see  [26],  [27],  [28],  [37],  [44]).  Moment  method  approximations  are  obtained 
from  truncated  Taylor  series  expansions  about  the  expected  value  of  the  input  parameter.  For  example, 
consider  a  function,  m(^),  expanded  about  the  mean  value,  C  The  first-order  accurate  approximation  for  the 
expected  value  of  u,  is: 


(2.9) 


1EfoN(0]  =  ^(0- 


Note  that  the  first-order  first  moment  (FOFM)  approximation  is  nothing  more  than  the  pointwise  (or  deter¬ 
ministic)  value  evaluated  at  the  mean  of  the  input.  Frequently,  this  is  referred  to  as  the  deterministic  solution. 
The  second-order  first  moment  (SOFM)  requires  the  computation  of  the  second  (sensitivity)  derivative  to 
improve  the  estimate  of  the  mean, 

(2.10)  TEsoHO]  =  uiO  +  ^Var(0  ^  • 

For  some  problems  we  investigated,  the  improvement  due  to  the  higher  order  correction  term  was  significant 
whereas  for  other  problems  it  was  not.  Estimates  of  the  variance  are  obtained  similarly.  The  first  order 
approximation  to  the  variance  of  u  is: 


(2.11) 


VarFo[w(6 


2 

Var(0. 


The  second  order  estimate  of  the  variance  of  u  is: 


(2.12) 


Varso[w(6 


2 

Var(0  +  \ 


For  discussion  purposes,  we  shall  refer  to  these  first-  and  second-order  approximations  to  the  second  moment 
as  FOSM,  SOSM  methods  respectively.  Extension  to  multiple  random  variables  is  straightforward  through 
Taylor  series  expansions  of  functions  involving  multiple  variables.  For  example,  ifu  =  u{^i,  ^2),  the  first-order 
moment  approximations  to  the  expected  value  and  variance  respectively  are 


(2.13) 


1Efo[w(^1,  ^2)]  = 


(2.14)  VarTO[w(^i,^2)]  =  (  ^ 


crl  + 


du 

9^2 


+  2 


du 


du 

9^2 


Covar(^i,^2)- 


where  the  covariance  (a  measure  to  the  extent  that  and  ^2  vary  jointly)  between  the  random  variables 
and  ^2  can  be  defined  in  terms  of  expected  values  as 


(2.15) 


Covar(^i,  ^2)  =  E  [^1^2  -  lE(^i)lE(^2)]  ■ 
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Table  2.2 

Quantiles  of  the  Standard  Normal  Distribution. 


a  (%) 

Quantile  (i^) 

67.000 

0.97 

90.000 

1.64 

95.000 

1.96 

99.000 

2.58 

99.900 

3.29 

99.990 

3.89 

99.999 

4.42 

Frequently,  it  is  useful  to  define  and  use  the  correlation  coefficient 

Covar(^i,^2) 


= 


1,11. 


The  extension  of  Eq.(2.14)  to  n  independent  random  variables  with  standard  deviations  J  results 

in  the  standard  error  estimate  for  (t„, 


(2.16) 


Note  the  similarity  of  this  FOSM  estimate  assuming  uncorrelated  data  with  the  propagation  of  error  formula, 
Eq.(2.3).  In  a  probabilistic  approach,  one  commonly  uses  an  input  uncertainty  with  a  specified  confidence 
interval,  e.g.  T  ±  1.96(T7jT,  a  95%  mean  confidence  interval.  Recall  that  a  confidence  interval  gives  a  bound 
within  which  a  parameter  is  expected  to  lie  with  a  certain  probability.  It  is  common  practice  to  distinguish 
between  single  prediction  confidence  intervals  in  which  a  randomly  drawn,  individual  (i.e.,  single)  sample  is 
expected  to  lie  with  a  pre-specified  probability  and  mean  confidence  intervals  in  which  the  population  mean 
is  expected  to  be  contained,  again  with  a  prescribed  probability.  Confidence  intervals  can  be  obtained  from 
knowledge  of  the  quantiles  of  the  distribution.  A  quantile  is  a  measure  of  the  location  of  a  point  such  that  a 
specified  fraction  of  the  data  lies  to  its  left.  For  example,  the  median  is  the  quantile  measuring  the  location 
of  the  point  such  that  50%  of  the  data  lie  to  its  left.  Quantiles  of  the  standard  normal  distribution  for  various 
confidence  intervals  are  shown  in  Table  2.2.  The  important  difference  between  the  interval  estimate  Eq.(2.3) 
used  by  Turgeon  and  the  probabilistic  approach,  Eq.(2.16),  is  the  interpretation  of  the  input  uncertainty 
and  hence  the  scale  factor  for  the  sensitivity  derivatives.  Note  that  ±3  standard  deviations  contains  more 
than  99%  of  the  interval  and  ±4  standard  deviations  contains  99.99%  of  the  interval.  If  experimental  data  is 
available,  then  estimates  of  the  input  standard  deviations,  . ,  can  be  obtained  that  will  be  less  (and  possibly 
substantially  less)  than  the  width  of  the  uncertainty  interval,  A^j,  used  in  Eq.(2.3).  This  distinction  is  further 
highlighted  in  Section  4  for  a  non-linear  convection-diffusion  problem.  Certainly,  repeated  experimental 
measurements  will  give  rise  to  a  probability  density  function  that  describes  the  statistics  of  the  distribution 
which  can  then  be  used  to  generate  accurate  and  reliable  probabilistic  output  for  multi-disciplinary  risk-based 
design  activities. 


2.5.  Polynomial  Chaos.  Recently,  several  papers  have  appeared  in  the  literature  investigating  spectral 
representations  of  uncertainty  (see  [14],  [15],  [16],  [17],  [47],  [48]).  An  important  concept  of  this  approach  is 
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the  decomposition  of  a  random  function  (or  variable)  into  separable  deterministic  and  stochastic  components. 
Specifically,  for  a  velocity  field  with  random  fluctuations,  we  write. 


(2.17) 


P 

i=o 


where  Ui{x)  is  the  deterministic  part  and  d>i(^)  is  the  random  basis  function  corresponding  to  the  mode. 
Effectively,  Ui{x)  is  the  amplitude  of  the  fluctuation.  The  discrete  sum  is  taken  over  the  number  of 
modes  represented,  P  =  which  is  a  function  of  the  order  of  the  polynomial  chaos,  p,  and  the  number 

of  random  dimensions,  n.  Here,  we  use  multi-dimensional  Hermite  polynomial  basis  functions  to  span  the 
n-dimensional  random  space  that  we  wish  to  represent  although  the  use  of  many  other  basis  functions  are 
possible.  A  convenient  form  of  the  Hermite  polynomials  is  given  by 


(2.18) 


—  ^  2  ^  ^ 


(-1)' 


dr, 


where  ^  =  (^ , . . . ,  ^ ^  J  is  the  n-dimensional  random  variable  vector.  As  discussed  in  [47] ,  there  is  a  one-to-one 
correspondence  between  the  functions  •••,  and  4>i(^).  The  Hermite  polynomials  form  a  complete 

orthogonal  set  of  basis  functions  in  the  random  space.  In  terms  of  the  inner  product, 

/OO 

fiOaiOwiOd^ 

-OO 


with  the  weight  function  W(^)  taking  the  form  of  an  n-dimensional  Gaussian  distribution  with  unit  variance 
defined  by 

(2.20)  ,rK)  ^ 

the  inner  product  of  the  basis  functions  is  zero  with  respect  to  each  other,  i.e.. 


(2.21) 


where  8ij  is  the  Kronecker  delta  function.  Once  the  modes,  Wfc,  of  the  solution  are  known,  then  statistics 
of  the  distribution  can  be  readily  evaluated.  The  mean  of  the  random  solution  is  given  by  ]Epc'[m]  =  uq. 
The  k  =  1,  ...,n  modes  are  the  Gaussian  estimates  of  the  variance,  all  higher  modes  provide  non-Gaussian 
interactions.  The  variance  of  the  distribution  is  given  by 

(2.22)  VarpcM  =  E  ■ 

i=l 


3.  Linear  Convection  with  a  Sonrce  Term.  The  primary  focus  of  this  example  is  on  the  application 
of  interval  analysis  for  both  time-dependent  and  steady-state  calculations.  We  consider  the  scalar  wave 
equation  with  a  source  term  designed  to  mimic  chemically-reacting  flow  problems  in  which  t  plays  the  role 
of  a  chemical  time  scale,  t  ^  1/^/,  where  kf  is  the  forward  reaction  rate.  The  governing  equation  is 


(3.1) 


du  du  u 

dt  ^  dx  T  ’ 


in  which  a  is  the  convective  wave  speed  (taken  to  be  a  =  1  for  all  cases  considered).  Near  chemical 
equilibrium,  reaction  rates  are  essentially  instantaneous,  i.e.,  kf  —>■  oo  and  t  — >  0,  resulting  in  disparate 
time  scales  that  cause  numerical  “stiffness”.  This  problem  has  been  studied  by  Godfrey  [18]  in  connection 
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Fig.  3.1.  Exact  time-dependent  solution  to  the  model  problem. 


with  implicit  preconditioning  algorithms.  Here,  we  consider  uncertainty  in  the  chemical  time  scale,  t,  due 
to  uncertainty  in  measured  reaction  rates.  The  exact  solution  to  this  equation  is 

(3.2)  u{x,t)  =  . 

a 

In  order  to  complete  the  definition  of  this  problem,  the  initial  and  boundary  condition,  respectively,  are 
taken  to  be 


(3.3) 

(3.4) 

With  these  conditions,  the  function,  g,  is 


u(x,  0)  =  1 
m(0,  t)  =  1 


X  e  [0,1], 
V  t  >  0. 


(3.5) 


git 


= 


1, 


t  >  i 

—  a 


o-{t-xla)lT  ^  0  <  t  <  -. 


Taking  the  chemical  time  scale,  r  =  0.2,  and  the  wave  speed,  a  =  1,  a  graph  of  the  exact  deterministic 
time-dependent  solution  from  t  =  0  to  t  =  1  is  shown  in  Figure  3.1.  For  the  uncertainty  analysis,  we  will 
examine  the  time-dependent  behavior  at  x  =  1  (seen  along  the  front  side)  and  the  steady  state  solution  (seen 
along  the  right  edge). 


3.1.  Interval  Analysis.  First-order  upwind  differences  were  used  to  approximate  the  spatial  derivative 
in  Eq.  (3.1).  The  outflow  boundary  condition  [x  =  1)  was  prescribed  by  approximating  =  0  to  first 

order  accuracy.  Three  time  integration  methods  were  implemented:  Euler  explicit,  Euler  implicit,  and  4- 
stage  Runge-Kutta.  In  terms  of  the  steady-state  residual,  Riu)  =  the  Euler  explicit  method  for 

ut  +  Riu)  =  0  can  be  written  as: 

(3.6)  -  AtR(w(”))  (Euler  Explicit-1). 


For  a  deterministic  problem,  that  is  all  there  is  to  this  method  (referred  to  hereafter  as  Euler  Explicit  -  1). 
However,  consider  the  case  in  which  there  is  uncertainty  in  the  input  chemical  time  scale,  t.  Let  t  be  defined 
to  be  the  interval 


(3.7) 


T  =  r[l-£,l-ke], 
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where  t  is  the  midpoint  of  the  interval  with  uncertainty  e.  The  residual  in  Eq.(3.6)  is  now  an  interval  rather 
than  a  point  value  at  any  x  location.  Moreover,  is  also  an  interval  after  the  first  iteration  as  a  consequence 
of  interval  operations  during  the  previous  time  step  even  if  the  initial  condition  has  no  uncertainty.  As  seen 
in  Figure  3.2,  this  method  results  in  exponential  growth  of  the  intervals  during  a  time  accurate  simulation. 
One  can  consider  alternatives,  namely. 


(3.8) 

«(”+!)  =  «(”)  -  AtR(il(”)) 

(Euler  Explicit-2) 

(3.9) 

«(”+!)  =«(")  -  AtR(M(”)) 

(Euler  Explicit-3) 

(3.10) 

«(”+!)  =«(")  -  AtR(il(”)) 

(Euler  Explicit-4) 

where  u  is  the  midpoint  of  the  interval.  Methods  2  and  4  evaluate  the  residual  with  the  midpoint  of  the 
interval  from  the  previous  time  step.  However,  is  still  an  interval  since  the  chemical  time  scale 

is  uncertain.  Since  Method  4  also  uses  the  midpoint  value  in  the  time  derivative,  it  gives  only  the  local 
uncertainty  at  each  step  due  strictly  to  the  uncertainty  in  t.  For  this  problem,  this  results  in  very  small 
interval  estimates  in  which  the  upper  and  lower  interval  bounds  visually  appear  to  be  on  top  of  one  another. 
Method  2  uses  the  interval  in  the  time  derivative  and  hence  accounts  for  cumulative  growth.  In  fact,  at 
any  point  in  time,  the  cumulative  (integrated)  sum  to  time  t  of  the  uncertainty  from  Method  4  is  equivalent 
to  the  Method  2  uncertainty  at  time  t.  Methods  1  and  3  both  use  the  interval  to  evaluate  the  residual 
and  hence  allow  for  temporal  growth  of  uncertainty.  The  magnitude  of  the  output  intervals  will  grow  without 
bound  in  Method  1  and  yet  with  Method  3  converge  to  steady-state  values  (provided  a  steady-state  exists). 
The  results  of  these  four  methods  with  81  equally-spaced  grid  points  along  the  x  —  axis  and  using  a  time 
step  At  corresponding  to  A  =  are  shown  in  Figure  3.2. 

It  is  common  practice  to  monitor  the  convergence  of  a  solution  process  in  CFD  simulations  in  terms  of 
a  vector  norm  of  the  steady-state  residual,  In  this  case,  however,  R{u^'^'>)  is  an  interval  not  suitable 

for  obtaining  point- wise  vector  norm  values  hence  we  monitored  the  L2  —  norm  of  R{u^'^'>).  The  residual 
histories  of  the  four  methods  are  shown  in  Figure  3.3.  For  Method  1,  the  interval  size  rapidly  reaches  the 
maximum  machine  representation,  on  this  computer  {  —  1.79769  x  10^°®,  1.79769  x  10®°®)  which  results  in  a 
midpoint  residual  of  zero  although  the  solution  is  completely  useless.  The  other  three  methods  smoothly 
converge  to  the  steady-state  in  t  <  2  yet  the  exact  steady-state  is  achieved  at  t  =  1.  However,  in  all  cases, 
the  residuals  at  t  =  1  have  been  reduced  more  than  three  orders  of  magnitude  and  the  values  at  this  time 
are  close  to  their  steady-state  values. 

One  of  the  most  common  4-stage  Runge-Kutta  methods  (and  the  one  implemented)  can  be  written  as 

7i=R(^i(”)), 

72  =  ^(«(”)+7i/2), 

(3.11)  73  =  i?(^;”;+72/2), 

74  =  ^(«(”)+73), 

Ar,(«)  =  _^(^^+272  +  273  +  74), 

«(”+!)  =«(")  -y 

Again,  we  face  the  same  issues  as  before,  namely  the  residual  evaluation,  R{u)  vs  R{u),  and  the  update  step, 
vs  Numerical  experiments  confirm  that  the  use  of  R{u)  for  the  residual  evaluations  and  in 

the  update  step  yields  the  most  useful  interval  analysis  results. 

For  many  CFD  simulations,  implicit  time  integration  methods  are  preferred.  For  demonstration  pur- 
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Fig.  3.2.  Interval  analysis  results  using  different  variations  on  the  Euler  Explicit  method. 
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Fig.  3.3.  Convergence  histories  for  the  four  implementations  of  the  Euler  Explicit  method  for  interval  analysis. 
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poses,  we  implemented  the  Euler  implicit  time  integration  scheme  which  in  delta  form  can  be  written  as 


(3.12) 


This  results  in  the  tri-diagonal  set  of  equations 


bi  Cl 

(  K  \ 

Au\ 

R{ur) 

a2  62  C2 

0 

Au2 

R{u2) 

(3.13) 

0 

1  ^n—1  ^n—l 

< 

AUn-l 

>  =  -  < 

\ 

dn 

AUn 

R(Un^ 

where  for  the  general  case  of  wave  speeds  of  either  sign, 


J 


At  the  inflow  boundary,  ci  =  R{ui)  =  0  and  6i  =  1  and  at  the  outflow,  an  =  — l.&n  =  ^,R{un)  =  0. 
For  positive  wave  speeds,  the  above  matrix  reduces  to  a  lower  bi-diagonal  matrix  which  can  be  solved 
by  a  forward  substitution  pass.  Likewise,  for  negative  wave  speeds,  the  matrix  (with  appropriate  change  in 
boundary  conditions)  reduces  to  an  upper  bi-diagonal  matrix  which  can  be  solved  by  a  backward  substitution 
step. 

Using  the  interval  midpoint  for  residual  evaluations,  a  comparison  of  the  time  history  from  the  Euler 
explicit,  Euler  implicit,  and  4-stage  Runge-Kutta  methods  at  a;  =  1  is  shown  in  Figure  3.4.  Although  difficult 
to  see,  the  Euler  implicit  and  Euler  explicit  methods  yield  identical  results  in  terms  of  interval  size.  The 
Runge-Kutta  method  results  in  larger  intervals  due  to  the  increased  number  of  residual  evaluations  and 
intermediate  update  steps  in  the  multi-step  method. 

If  one  is  only  interested  in  the  steady-state  solution,  other  options  become  available,  primarily  through 
the  use  of  algorithms  that  are  stable  for  large  time  steps  and  in  some  cases,  space-marching  methods.  In  this 
case,  both  are  applicable.  For  example,  the  Euler  implicit  method,  Eq.  (3.12)  can  be  used  with  At  oo. 
Since  the  governing  equation  is  linear,  the  steady  solution  is  obtained  in  one  step  by  solving  the  resulting 
linear  problem  which  we  refer  to  as  a  Time  Alarch  method.  The  other  alternative  is  to  implement  a  Space 
Alarch  method  by  eliminating  the  time  derivative  in  Eq.  (3.1),  i.e.. 


(3.15) 


u 

T 


Discretizing  this  equation  with  first-order  upwind  differences  (taking  a  >  0)  yields 


(3.16) 


where  Uj  =  u{jAx).  The  results  of  the  Time  Alarch  versus  Space  Alarch  approaches  can  be  seen  in  Figure 
3.5  where  the  advantage  of  space  marching  is  evident. 


12 


Fig.  3.4.  Twie  accxirate  interval  calculations  with  Euler  explicit,  Runqe-Kutta,  and  Euler  implicit  time  integration  methods. 


Fig.  3.5.  Steady  state  interval  sohUion  usmg  the  Euler  implicit  method  with  an  infinite  time  step  (Time  March)  cornpared 
to  a  space-marching  method. 
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Fig.  3.6.  Histograms  of  the  mput  chemical  time  scale,  r,  and  the  output  variable,  u  at  x  =  1  from  1000  Monte  Carlo 
simulations. 
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Fig.  3.7.  Mean  values  and  95%  single  prediction  confidence  intervals  from  time-dependent  Monte  Carlo  at  three  x  locations. 


Fig.  3.8.  Comparison  of  the  First-Order  Moment  method  with  1000  Monte  Carlo  simulations. 


3.2.  Monte  Carlo  and  Moment  Method.  The  Monte  Carlo  and  First-Order  moment  methods  were 
also  implemented.  For  these  probabilistic  methods,  we  took  the  chemical  time  scale,  t,  to  be  a  Gaussian 
random  variable  with  a  mean,  t  =  0.2  and  a  5%  coefficient  of  variation  (CoV).  Figure  3.6  shows  histograms 
from  1000  Monte  Carlo  simulations  where  it  can  be  seen  that  the  output  is  clearly  non-Gaussian. 

Time-dependent  Monte  Carlo  calculations  are  compared  at  three  x  locations  in  Figure  3.7.  For  this 
problem,  the  solution  at  any  x  decays  from  the  initial  condition  of  w  =  1  until  its  steady-state  value  is 
reached  at  t  The  95%  mean  confidence  intervals  are  roughly  30  times  narrower  and  are  not  shown 

here  for  clarity.  Steady-state  results  from  the  FO  moment  method  are  compared  with  1000  Monte  Carlo 
simulations  in  Figure  3.8.  The  sensitivity  derivative  required  to  implement  the  FO  moment  method  is 


(3.17) 


^  ^  ^  -x/ar 

dr 


The  95%  confidence  intervals  from  the  FOSM  are  comparable  to  the  single  prediction  95%  confidence  intervals 
from  Monte  Carlo.  The  distribution  of  standard  deviation  and  coefficient  of  variation  are  in  good  overall 
agreement. 
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Fig.  4.1.  Exact  deterministic  solution  and  first  three  sensitivity  derivatives. 


4.  Non-linear  Burgers  Equation.  Computational  fluid  dynamicists  frequently  study  simplified  prob¬ 
lems  designed  to  mimic  certain  features  of  a  more  complicated  situation  but  at  lower  cost  and  effort.  Many 
such  problems  exist  but  one  that  is  particularly  useful  is  the  model  non-linear  convection-diffusion  problem 
proposed  by  Rakich  and  referred  to  as  the  general  Burgers  equation  [4].  Based  on  specific  choices  of  para¬ 
meters  in  the  equation,  one  can  obtain  linear  or  non-linear  behavior.  We  consider  the  steady  form  of  this 
equation,  which  in  conservation  law  form  is 


(4.1) 


df  d  f  du\ 
dx  dx  \  dx ) 


R{u)  =  0 


where  the  flux,  /,  is  a  non-linear  function  of  u.  We  took  the  specific  case  of  /  =  w  (^  —  w)  for  which  the 
exact  stationary  solution  is  given  by 

(4.2)  u{x)  =  ^  l-ktanh(^^^  . 

For  the  uncertainty  analysis,  we  took  the  input  viscosity  to  be  stochastic,  hence  the  sensitivity  derivatives, 
d^u/dijX,  are  important  in  the  analysis.  The  exact  solution,  u{x),  and  its  first  three  derivatives  with  respect 
to  viscosity,  fi,  are  shown  in  Figure  (4.1).  The  solution  rises  monotonically  from  0  at  — cx)  to  1  at  -|-cx).  Note 
that  the  higher  derivatives  are  increasingly  oscillatory.  One  practical  consequence  of  this  that  we  found  was 
the  need  for  high-order  methods  or  finer  grids  to  numerically  estimate  those  derivatives. 

For  all  numerical  computations,  we  limited  the  computational  domain  to  x  E  [—3,  3].  Truncating  infinite 
domains  in  CFD  is  commonplace,  external  flow  over  airfoils,  wings,  and  aircraft  configurations  are  a  few 
examples.  However,  there  are  consequences  of  this  with  regard  to  boundary  condition  treatment  of  stochastic 
variables.  A  more  complete  discussion  of  this  topic  and  a  random  held  analysis  of  this  problem  can  be  found 
in  [29]. 
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Fig.  4.2.  Numerical  solutions  using  a  second-order  accurate  centered  scheme  on  three  mesh  levels  (left).  Typical  conver¬ 
gence  history  using  Newton’s  method  for  the  general  Burger  equation  (right). 


4.1.  Deterministic  Problem.  Newton’s  method  was  used  to  solve  the  non-linear  problem,  Eq.(4.1). 
This  results  in  the  linearized  system  of  equations  and  the  update  step 


(4.3) 


«(”+!)  =«(")  + 


Second-order  centered  differences  were  used  to  approximate  the  spatial  derivatives.  Dirichlet  boundary 
conditions  were  specified  at  the  endpoints,  x  =  ±3,  from  the  exact  deterministic  solution.  This  results  in  a 
tri-diagonal  set  of  equations  to  solve  for  each  Newton  iteration  (for  details,  see  [29] ) 


bi  Cl 

f  A  > 

Awi 

R{ui) 

<32  b2  C2 

0 

Au2 

R{u2) 

(4.4) 
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where 


(4.5) 


(  n  -  _ ^ 

“  2Ax  Aa;2 


b  ■  — 

.  _  (I/2-M3-  +  1) 
'-'3  ~  2Ax 


bi  =  bn  =  1,  and  ai  =  ci  =  a„  =  c„  =  R{ui)  =  R{un)  =  0. 

Numerical  solutions  were  generated  on  a  sequence  of  uniformly-spaced  grids  for  which  the  number  of 
grid  points  and  mesh  spacing  is  given  in  Table  4.1.  The  numerical  solution  for  the  dependent  variable,  u,  on 
three  different  grids  is  shown  in  the  left  of  Figure  4.2.  A  typical  convergence  history,  corresponding  to  the 
129  mesh  point  solution  is  on  the  right.  Independent  of  the  level  of  grid  refinement,  the  residual  converged 
to  machine  epsilon  in  the  same  number  of  iterations  (±1). 

Several  approaches  have  been  proposed  for  establishing  truncation  error  estimates.  Probably  the  most 
common  approach  is  based  on  Richardson  extrapolation  although  other  estimates,  such  as  mult-grid  estimates 
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Fig.  4.3.  Left  -  Convergence  of  the  L2  error  norm  versus  mesh  spacing.  Bight  -  Distribution  of  exact  and  estimated 
discertization  error  on  three  grids. 


can  also  be  used.  Roache  [40]  has  used  Richardson  extrapolation  with  modifications  to  obtain  an  error 
estimator  and  this  technique  will  be  used  here  for  demonstration  purposes.  The  discretization  error  estimates 
that  we  require  can  be  obtained  by  series  expansions.  The  result  is: 

(4-6) 

where 


E1.  E2  =  error  estimates  on  the  fine  grid  “1”  and  coarse  grid  “2”,  respectively 

£  =  U2  —  ui;  the  pointwise  difference  between  successive  solutions 

r  =  —  >  1;  the  ratio  of  mesh  spacing 
hi 

p  =  the  order  of  accuracy  of  the  numerical  method. 

The  Grid  Convergence  Index  (GCI)  used  by  Roache  as  a  discretization  error  estimator  is  defined  by 
(4.7)  GCI  =  Es\Ei\  or  GCI  =  Es\E2\ 

where  Es  is  a  factor-of-safety  with  a  recommended  value  of  Es  =  1.25  (which  was  used  here)  for  high  fidelity 
numerical  experiments.  In  all  cases  examined,  the  GGI  estimate  bounded  the  true  discretization  error  at 
every  point  in  the  domain.  Note  that  error  estimates  are  available  on  all  grid  levels.  Roy  [43]  has  extended 
this  analysis  to  mixed-order  schemes  and  further  refinements  are  possible  using  information  from  additional 
grids. 

Measures  of  accuracy  and  discretization  error  are  shown  in  Figure  4.3.  The  plot  on  the  left  shows 
L2  —  norm  of  the  actual  error.  The  slope  of  the  line  is  well  known  to  be  a  measure  of  the  order  of  accuracy 
of  the  discretization.  Table  4.1  shows  the  computed  values  of  the  slope  of  each  line  segment.  The  exact 
distribution  of  the  discretization  error  is  shown  on  the  right  of  Figure  4.3  along  with  the  Grid  Convergence 
Index  (GCI)  estimate. 
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Table  4.1 

Summary  of  grid  convergence  parameters  and  results  for  Burger’s  eguation. 


Grid  Level 

Number  of 

Grid  Points 

Mesh 

Spacing,  Ax 

L2  norm 

of  Error 

Slope  of 
Segment 

1 

9 

0.75 

0.01383131 

- 

2 

17 

0.375 

0.00350755 

1.97940 

3 

33 

0.1875 

0.00095693 

1.87398 

4 

65 

0.09375 

0.00023901 

2.00129 

5 

129 

0.046875 

0.00006023 

1.98851 

6 

257 

0.0234375 

0.00001505 

2.00001 

7 

513 

0.0117188 

0.00000376 

2.00000 

4.2.  Stochastic  Problem. 


4.2.1.  Exact  Stochastic  Solution.  The  viscosity,  ft,  was  taken  to  be  a  Gaussian  random  variable, 
PAiild),  with  a  mean  Jl  =  0.25  and  a  coefficient  of  variation,  CoVift)  =  —  =  10%.  The  exact  stochastic 
solution  can  be  expressed  in  terms  of  its  probability  density  function,  which  is  given  by 

-1 


(4.8) 


pu{u{x))  = 


dpi 


u{x,ij) 


PMild). 


Justification  for  this  expression  is  rather  intuitive  since  it  is  simply  a  statement  that  the  probability  of  a 
random  variable  following  within  a  specified  range  is  constant  under  a  transformation  of  variables.  Using 
Eq.(4.1)  and  the  definition  of  the  Gaussian  PDF,  this  expression  can  be  analytically  evaluated  to  obtain 


•  (l-2u 


(4.9) 


pu{u{x))  = 


(li— l)iitanh  ^  (1  — 


8\/^(t 

A  plot  of  this  function  is  shown  in  Figure  4.4  followed  by  the  enlarged  view  in  Figure  4.5  with  more  resolution 
of  the  PDF.  It  can  be  seen  that  pu  varies  markedly  with  x.  Near  the  boundaries,  the  PDF  is  highly  skewed 
and  peaked  but  approaches  a  Gaussian  as  a;  ^  ±1  from  the  nearest  boundary.  At  a;  =  0,  the  solution  is 
deterministic,  i.e.,  the  variance  is  zero,  which  all  random  variable  models  predict.  As  a  side  note,  random 
field  models  may  produce  quite  different  results,  particularly  near  the  centerline  [29]. 

4.2.2.  Interval  Analysis.  We  applied  interval  analysis  to  both  the  exact  solution  and  the  GFD- 
oriented  numerical  method  described  previously.  The  input  interval  object  was  the  viscosity,  which  was 
defined  to  be 


(4.10)  pL  =  /l[1  —£,!  +  £]. 

Numerical  results  obtained  by  evaluating  the  exact  deterministic  solution,  Eq.  (4.2)  for  an  input  value, 
£  =  0.1,  are  shown  in  Figure  4.6.  The  error  bars  indicate  the  width  of  the  output  interval.  One  of  the 
most  attractive  things  about  interval  analysis  is  that  it  is  very  simple  to  implement.  However,  as  seen  in  the 
interval  analysis  results  from  the  numerical  method  in  which  both  the  Jacobian  matrix  and  residual  were 
evaluated  with  the  interval  midpoint  values.  Figure  4.7,  the  uncertainty  estimates  are  unacceptably  large. 
The  bounds  became  so  large  for  the  bottom  two  cases  that  a  change  in  scale  was  required  to  keep  the  error 
bars  within  the  display  area. 
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Fig.  4.4.  Probability  density  function  of  the  exact  stochastic  solution  of  the  Burger  equation. 


Fig.  4.5.  Enlarged  view  of  the  exact  probability  density  function  of  u  for  the  Burger  equation. 


Fig.  4.6.  Interval  analysis  of  the  exact  deterministic  solution  for  an  input  uncertainty  in  the  viscosity,  e  =  0.1(2.e.l0%) . 


19 


Fig.  4.7.  Interval  analysis  results  from  the  numerical  solution  of  Burger’s  equation.  Note  the  change  in  scale  on  the 
bottom  two  plots. 

4.2.3.  Moment  Methods.  The  moment  methods  require  the  sensitivity  derivatives  with  respect  to 
the  input  random  variable.  In  this  example,  we  computed  numerical  approximations  from  the  continuous 
sensitivity  equation  (CSE)  approach  and  the  discrete  adjoint  (DA)  method.  A  brief  description  of  these  two 
techniques  follows: 

Continuous  Sensitivity  Equation  (CSE).  The  conservation  law  form  of  the  Burgers  equation  (4.1)  can 
be  re-written  in  quasi-linear  form  as 

(4.11)  Xux  =  lauxx  where  X  =  —  —  u. 

The  continuous  sensitivity  equation  is  derived  by  differentiating  Eq.  (4.11)  with  respect  to  noting  that 
u  =  u{x;ij,).  Defining  =  du/dji  and  carrying  out  the  differentiation  yields  a  linear  differential  equation 
for  the  sensitivity  derivative,  namely, 

(4.12)  ~  T  Uxx- 

Given  a  numerical  solution  of  the  Burgers  equation  as  input,  one  can  use  a  variety  of  numerical  methods  to 
solve  Eq.  (4.12).  We  used  a  spatial  discretization  consistent  with  the  flow  solution  method  (i.e.,  a  2"^*  order 
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accurate  centered  scheme)  and  Newton’s  method.  Since  the  problem  is  linear,  one  Newton  iteration  always 
finds  the  solution  to  machine  precision  («  10^^®).  It  is  worth  mentioning  that  the  Jacobian  matrix  of  the 
CSE  equation  is  identical  to  the  Jacobian  matrix  of  the  Burgers  equation.  This  simplifies  the  coding  of  the 
CSE  since  one  merely  has  to  change  the  right  hand  side  in  the  original  software. 

Discrete  Adjoint  Method.  The  discrete  adjoint  method  is  commonly  used  in  aerodynamic  optimization 
problems  ([34])  in  which  a  user  defined  cost  function  is  minimized.  The  cost  function,  denoted  F(u]  ?i)  where 
represents  the  generic  design  variable,  is  augmented  with  a  vector  of  Lagrange  multipliers  operating 
on  the  discrete  version  of  the  governing  equations,  in  this  case  Eq.  (4.1),  i.e., 

(4.13)  F*  =  F{u;^i)  +  \^R{u;u). 


This  equation  is  next  differentiated  with  respect  to  the  design  variable,  or  more  generally,  any  generic  variable 
denoted  After  some  rearrangement,  one  obtains 


(4.14) 


dF*  _  ^  xT  (^\ 


Since  the  vector  of  Lagrange  multipliers  is  arbitrary,  the  first  term  in  parentheses  may  be  equated  to  zero  to 
yield  a  linear  system  of  equations  for  A, 


(4.15) 


f^Yx--— 


Once  the  vector  of  Lagrange  multipliers  is  known,  the  sensitivity  of  the  cost  function  with  respect  to  can 
be  evaluated  from 


(4.16) 


This  basic  implementation  is  widely  used  in  a  design  environment  (see  [34])  and  precludes  the  need  to 
compute  the  flow  sensitivities,  i.e,.  9^/9?^,  directly.  In  this  problem,  we  require  those  derivatives  and  they 
can  be  readily  obtained  by  defining  the  set  of  cost  functions,  F  =  {Fj  =  Uj\j  =  1,2,  ...n]  where  n  is  the 
total  number  of  grid  points.  Extending  Eq.  (4.16)  to  a  system  of  functions  and  inserting  the  definition  of  F 
yields 


(4.17) 


I 


where  /  is  the  identity  matrix  and  A  is  a  matrix  consisting  of  a  set  of  column  vectors,  \j  for  each  j.  The 
required  sensitivity  derivatives  can  now  be  found  from  Eq.  (4.14),  which  upon  substitution  and  rearrangement 
reduces  to 


(4.18) 


du  J,  f  dR\ 

(itej 


Note  that  Eq.  (4.17)  can  be  solved  efficiently  for  all  right  hand  sides  by  LU  decomposition.  All  sensitivity 
derivatives  can  then  be  computed  by  the  single  matrix- vector  multiply  indicated  in  Eq.  (4.18). 

Although  it  may  not  be  readily  apparent,  both  the  CSE  and  the  DA  methods  yield  identical  results 
provided  both  methods  are  formulated  and  implemented  consistently.  In  practice,  this  generally  does  not 
occur  because  of  differences  in  coding,  convergence  level,  etc.  Here,  however,  the  equations  are  always  solved 
to  machine  precision,  the  discretizations  are  consistent,  and  consequently,  the  sensitivity  derivatives  obtained 
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Fig.  4.8.  Left  -  First  sensitivity  derivative  computed  by  the  discrete  adjoint  and  continuous  sensitivity  methods  compared 
to  the  exact  solution.  Bight  -  Comparison  of  probabilistic  (dashed  lines)  vs.  non-pro babilistic  (error  bars)  interpretation  of 
input  uncertainty  on  x  G  [—3,0]. 


by  both  methods  are  identical.  A  plot  of  the  first  derivative  with  respect  to  viscosity  by  both  methods  is 
shown  on  the  left  of  Figure  4.8  and  compared  with  the  exact  derivative.  The  right  half  of  the  figure  shows 
results  for  estimating  the  uncertainty  assuming  that  1)  the  input  uncertainty  contains  100%  of  the  error 
(non-probabilistic  -  Eq.  (2.3))  and  2)  a  probabilistic  FOSM  approximation  assuming  that  the  input  error 
represents  three  standard  deviations  («  99%  of  the  samples).  The  advantage  of  knowledge  about  the  input 
data  is  clearly  seen  in  the  reduced  uncertainty  in  the  output. 

A  comparison  of  the  absolute  error  in  the  mean  (jl  =  0.25)  and  relative  error  in  the  standard  deviation 
{CoV {fx)  =  ^  =  10%)  predicted  by  Monte  Carlo  and  the  first-  and  second-moment  methods  are  shown 
in  Figure  4.9.  The  distribution  of  the  mean  error  from  the  first-order  moment  method  closely  follows  the 
second  derivative  distribution  shown  in  Figure  4.1.  Adding  the  second  derivative  correction  term  to  obtain 
the  second  order  moment  estimate  of  the  mean  results  in  substantially  less  error.  As  can  be  seen  on  the 
right,  the  relative  error  in  estimating  the  standard  deviation  is  approximately  6%  and  2%  for  the  FOSM  and 
SOSM  estimates  respectively  at  the  boundaries.  One  thousand  MC  simulations  results  in  roughly  the  same 
maximum  error  as  the  SOSM  method  at  roughly  300  times  the  cost.  Increasing  the  number  of  simulations 
an  order  of  magnitude  reduced  the  error  by  approximately  three  at  considerable  additional  cost. 


4.2.4.  Polynomial  Chaos.  Substituting  in  polynomial  chaos  expansions  for  the  dependent  variable, 
u,  and  the  input  random  variable,  fi,  into  the  residual,  Eq.(4.1),  and  taking  the  inner  product,  (•,  yields 
an  equation  for  the  residual  of  the  mode. 


(4.19) 


Rk 


ld_ 

2fe 
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i=0 j=0 


P  P 

EE 

i=0 


dx^ 


where  eijk  =  We  iteratively  solved  this  equation  by  the  Euler  explicit  method  but  found  that 

the  CPU  time  to  achieve  the  machine  zero  steady-state  solution  was  excessive  relative  to  implicit  methods. 
Consequently,  Newton’s  method  was  implemented  by  using  the  linearization  of  the  residual,  which  can  be 
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Fig.  4.9.  Absolute  error  in  the  mean  (left)  and  relative  error  in  the  standard  deviation  (right)  using  Monte  Carlo  and 
first-  and  second-moment  methods. 


shown  to  reduce  to 


(4.20)  Ai?fc 


/  p  p  \  ■ 

\  i=^  i=o  ) 


p  p 

-EE  ^ijkfdi 

i=0 


9^  Am; 
dx^ 


where  8ij  is  the  Kronecker  delta  function.  The  boundary  conditions  were  specified  from  the  exact  stochastic 
solution  by  matching  moments  of  the  distribution  up  to  a  user-specified  order.  The  three  boundary  conditions 
studied  were  obtained  by  matching 

1.  the  mean  (BCl) 

2.  the  mean  and  variance  (BC2) 

3.  the  mean,  variance,  and  skewness  (BC3) 

The  first  boundary  condition,  BCl,  implies  that  the  first  mode,  uq  is  set  to  the  mean  of  the  exact  solution 
and  the  remaining  modes  are  set  to  zero  on  the  boundary.  The  second  boundary  condition  sets  uq  and  mi 
to  the  mean  and  variance  of  the  exact  solution,  respectively,  the  remaining  modes  are  set  to  zero.  Similarly, 
BC3  sets  mq  to  the  mean  and  solves  a  simple  algebraic  problem  to  find  the  values  of  mi  and  U2  such  that  the 
variance  and  skewness  of  the  polynomial  chaos  distribution  match  the  exact  solution. 

One  advantage  of  Polynomial  Chaos  is  that  the  output  PDF’s  which  depend  on  the  order  of  the  chaos 
can  be  easily  obtained.  For  a  single  random  variable,  the  first-order  chaos  yields  a  Gaussian  distribution 
with  mean  uq  and  standard  deviation  |mi|,  i.e., 

-{u-Vfif- I2u\ 

(4-21)  VpcM)  =  - ^ — • 

V27r  \ui 

The  higher-order  modes  account  for  non-Gaussian  interactions  and  are  reflected  in  the  output  PDF.  For  a 
second-order  chaos,  the  result  is 


PPC2  {U)  = 


-  +  u\+iU2  (M-M0-I-M2))  /8m2 

/27r  \ul  +  4m2  (m  —  Mo  +  M2) 


Further  analytic  representations  for  higher  order  chaoses  are  possible  but  lengthy. 
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Fig.  4.10.  The  first  four  modes  from  Polynomial  Chaos.  Compare  to  the  sensitivity  derivatives  of  the  exact  solution. 


Numerous  results  were  generated  using  Hermite  polynomial  chaoses  of  varying  order  with  different 
boundary  condition  treatment  and  on  different  grids.  Here  we  summarize  these  results.  The  first  four  modes 
from  an  order  3  PC  are  shown  in  Figure  4.10.  Recall  that  mode  0  represents  the  expected  value  (i.e.,  the 
mean).  Note  the  similarity  of  the  shapes  of  the  higher  modes  with  the  sensitivity  derivatives  of  the  exact 
solution  in  Figure  4.1.  The  only  discrepancy  that  can  be  seen  is  between  the  third  mode  and  the  third 
derivative  which  is  due  to  boundary  condition  treatment  for  this  mode  and  will  be  discussed  subsequently. 

Figure  4.11  shows  absolute  errors  in  the  mean  and  relative  errors  in  the  standard  deviation  using  order 
1  and  2  PC  on  129  and  513  point  grids.  Refining  the  grid  for  the  first-order  chaos  gives  a  slightly  larger 
maximum  error  than  the  coarser  grid.  Computations  show  that  the  first-order  method  is  grid  converged 
on  the  129  point  mesh,  i.e.,  further  grid  refinement  makes  no  appreciable  reduction  in  error  nor  is  there 
significant  increases,  just  re-distribution.  However,  the  second-order  polynomial  chaos  is  clearly  not  grid 
converged  on  the  129  point  grid.  There  is  a  substantial  reduction  in  the  mean  error  on  the  finer  513  point 
grid.  In  terms  of  the  estimates  of  the  standard  deviation,  the  first-order  chaos  results  on  both  grids  are 
similar  and  very  close  to  the  SOSM  method  (see  Fig.  4.9).  The  2"^*  order  chaos  has  very  small  error  in  the 
standard  deviation,  almost  an  order-of-magnitude  smaller  than  the  estimate  obtained  from  10,000  Monte 
Carlo  simulations. 

A  further  examination  of  the  error  convergence  versus  the  order  of  the  PC  was  conducted.  In  theory, 
exponential  reduction  rate  of  the  error  as  the  order  of  the  chaos  is  increased  should  be  observed.  In  practical 
CFD  simulations,  there  are  issues  that  are  likely  to  make  such  convergence  rates  unachievable  although  this 
does  not  preclude  the  use  or  utility  of  the  polynomial  chaos  method.  Figure  4.12  supports  this  statement 
where,  on  the  left,  the  convergence  behavior  of  the  L2  norm  shows  that  in  order  to  achieve  the  “theoretical” 
convergence,  more  and  more  information  needs  to  be  supplied  at  the  boundary  as  the  order  of  the  chaos 
increases  (see  Table  4.2).  This  simply  means  that  in  order  to  match  the  exact  solution  to  a  given  level  of 
accuracy,  boundary  conditions  must  also  be  specified  consistent  with  the  level  of  accuracy  sought.  Thus,  to 
stay  on  the  exponential  convergence  curve,  higher  order  chaoses  require  higher  order  statistics  at  the  bound¬ 
aries  of  the  computational  domain.  Note  also  that  finer  grids  are  also  required  to  stay  on  the  exponential 
reduction  rate  curve  due  to  the  increased  frequency  content  inherent  in  the  higher  order  modal  solutions. 
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Fig.  4.11.  Absolute  error  in  the  mean  (left)  and  relative  error  in  the  standard  deviation  (right)  from  Polynomial  Chaos. 


Provided  that  the  PDF  of  the  output  uncertainty  on  the  boundary  is  known,  one  can  determine  the 
required  modal  values,  on  the  boundary  by  matching  moments  as  was  done  here.  It  is  not  likely  that  the 
PDF  will  be  available  for  general  CFD  simulations.  Moreover,  we  will  be  fortunate  to  get  good  estimates  of 
the  mean  and  variance  at  the  boundaries.  In  the  absence  of  higher  order  statistics,  one  is  forced  to  make 
assumptions  that  degrade  the  error  reduction  rate  of  the  method.  However,  as  shown  earlier,  even  relatively 
low-order  (<3)  chaos  can  produce  results  superior  to  the  first-  and  second-moment  methods. 

To  further  demonstrate  the  boundary  condition  issue,  the  right  of  Figure  4.12  shows  the  error  convergence 
for  a  slightly  simpler  problem,  the  1—D  heat  equation  with  internal  heat  generation.  The  significant  difference 
between  these  two  problems  is  the  stochastic  boundary  conditions.  For  the  heat  equation,  the  domain  is 
finite  and  the  boundary  conditions,  by  problem  definition,  have  no  uncertainty.  The  uncertainty  enters 
through  the  thermal  conductivity  which  was  treated  as  a  Gaussian  random  variable  with  a  CoV  =  10%. 
This  problem  has  a  simple  exact  stochastic  solution  and  since  there  is  no  boundary  uncertainty,  the  exact 
boundary  condition  treatment  is  to  set  all  higher  modes  to  zero.  As  seen,  the  exponential  convergence  rate 
is  achieved.  On  the  other  hand,  for  the  Burgers  equation,  the  effect  of  not  being  able  or  willing  to  specify 
the  higher  modes  exactly  results  in  the  “peeling  off”  of  the  error  from  the  theoretical  curve. 

Further  impact  of  the  boundary  condition  specification  for  the  Burgers  equation  can  be  seen  in  Figure  4.13 
in  which  the  exact  solution  is  compared  with  the  probability  density  function  of  a  third-order  chaos  at  a;  =  —2; 
the  shaded  areas  highlight  their  difference.  The  three  boundary  conditions  correspond  to  specification  of  the 
mean  (BCl),  mean  and  variance  (BC2),  mean,  variance,  and  skewness  (BC3)  respectively,  i.e.,  the  boundary 
condition  specifications  are  increasingly  complete.  The  first  boundary  condition  shows  that  ignoring  the 
variance  at  the  boundary  results  in  a  much  more  peaked  distribution  than  should  be  predicted.  This  effect 
is  even  more  pronounced  closer  to  the  boundary  on  which  there  is  no  distribution,  just  a  mean  value  and 
zero  variance.  Consequently,  moving  slightly  inward  results  in  some  but  little  variance  since  it’s  forced  to 
zero  nearby,  and  hence  a  large  narrow  peak.  However,  as  seen  in  the  lower  left  of  the  figure,  addition  of  the 
variance  makes  substantial  improvement  in  the  prediction  although  the  distribution  is  Gaussian.  Providing 
additional  information  concerning  the  boundary  skewness  is  seen  to  shift  the  distribution  markedly  such  that 
the  interior  skewness  is  now  well  predicted. 
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Table  4.2 

Case  description  for  Figure  4-12 


Case 

Boundary  Condition 

Grid  Points 

1 

Mean  only 

129 

2 

Mean  and  Variance 

129 

3 

Mean,  Variance,  and  Skewness 

129 

4 

Mean,  Variance,  and  Skewness 

257 

5 

Mean,  Variance,  and  Skewness 

513 

Fig.  4.12.  Error  in  the  mean  for  various  orders  of  Polynomial  Chaos.  Left  -  Burger’s  eguation  with  various  boundary 
conditions  and  grids.  Bight  -  Heat  eguation  with  stochastic  internal  heat  generation.. 


The  statistics  of  the  distribution  predicted  by  a  4*^  order  chaos  are  compared  with  the  exact  solution 
in  Figure  4.14.  Without  symbols  on  the  plots,  the  distributions  are  all  but  indistinguishable.  Note  that  the 
standard  deviation  is  a  maximum  in  the  vicinity  of  a;  =  ±1  where  the  PDF  is  nearly  Gaussian  as  indicated 
by  its  skewness  near  0  and  kurtosis  value  near  3.  The  PDF’s  switch  from  skewed  to  the  left  to  skewed  to 
the  right  in  each  half  domain  x  G  [—3,  0}  and  x  <E  {—3,  0].  It  is  also  clear  that  the  departure  from  Gaussian 
is  maximal  at  the  boundary. 

An  estimate  of  the  computational  effort  associated  with  this  problem  is  shown  in  Figure  4.15.  The 
deterministic,  first-  and  second-order  moment  methods,  and  polynomial  chaos  of  orders  1-3  are  compared 
with  the  cost  of  100  Monte  Garlo  simulations.  Moment  methods  cost  2-3  times  more  than  a  deterministic 
solution.  The  work  of  polynomial  chaos  is  relatively  high,  about  20  times  the  cost  of  a  deterministic  solution 
but  there  are  opportunities  to  reduce  this  cost  (e.g.  loosely-coupled  algorithms,  multi-grid,  ...). 
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Fig.  4.13.  Comparison  of  exact  and  polynomial  chaos  pdf’s  at  x  —  —2  with  three  different  boundary  conditions.  The 
shaded  areas  highlight  the  difference  between  the  exact  and  approximate  PDF’s. 


Fig.  4.14.  Statistics  of  the  output  distribution  from  4*^  order  Polynomial  Chaos  compared  with  the  exact  solution. 


Relative  CPU  Times 


Fig.  4.15.  Pelative  cpu  times  to  generate  solutions  to  the  Burger  eguation. 
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Kurtosis 


[3 ,  Shock  Wave  Angle 


Fig.  5.1.  Sketch  of  supersonic,  inviscid  flow  over  a  wedge. 


5.  Oblique  Shock  Waves.  Shock  waves  and  expansion  fans  are  fundamental  building  blocks  for  in¬ 
viscid  compressible  flow  theory.  We  begin  by  addressing  inviscid,  supersonic  flow  over  a  wedge  at  Mach 
numbers  and  wedge  angles  for  which  an  attached,  steady,  oblique  shock  forms.  A  sketch  of  the  problem  is 
shown  in  Figure  5.1.  The  output  quantity  of  interest  for  this  example  is  the  pressure  rise  across  the  shock 
wave,  which  is  strictly  a  function  of  the  upstream  Mach  number.  Mi  and  the  wedge  angle  9  (for  a  fixed 
7).  In  all  cases,  we  take  the  ratio  of  specific  heats,  7  =  1.4. 

5.1.  Deterministic  Problem.  For  a  perfect  gas,  the  pressure  rise  across  a  normal  shock  wave  is  given 
by  the  Rankine-Hugoniot  relation 

(5.1)  ^  =  1  +  ^(Mf  -  1) 

Pi  7+1 

where  the  subscripts “1”  and  “2”  refer  to  conditions  just  ahead  of  and  behind  the  shock  wave  respectively. 

The  pressure  rise  across  an  oblique  shock  can  be  obtained  from  the  normal  shock  relation  by  replacing 
the  Mach  number  in  Eq.(5.1)  by  the  normal  component  of  the  upstream  Mach  number,  Mi,^  in  this  case. 
From  geometric  considerations. 


(5.2) 


Mi,^  =  Ml  sin  /3. 


A  relationship  for  the  shock  wave  angle,  /3,  can  be  obtained  from  the  tangential  momentum  equation  (the 
so-called  (3  —  9  —  M  relationship).  A  common  form  of  this  relationship  is 


(5.3) 


tan  9  =  2  cot  j3 


m2  sin2  (3-1  \ 

M2  (7 -y  cos  2/3)  +  2  J 


For  a  detailed  derivation,  see  e.g.  [5].  Note  that  this  relationship  uniquely  specifies  the  flow  deflection 
angle,  9,  for  a  given  Mach  number  ahead  of  the  shock,  (denoted  generically  by  M)  and  shock  wave  angle,  (3. 
However,  in  practice,  one  generally  knows  the  Alach  number  and  the  geometry,  9,  (to  within  a  manufacturing 
tolerance),  thus,  one  needs  to  solve  Eq.(5.3)  for  (3.  With  some  manipulation,  this  equation  can  be  written 
as  a  cubic  polynomial  in  the  variable  sin2  (3,  hence  three  solutions  for  (3  exist  for  each  {M,  9}  pair.  Two  of 
the  solutions  correspond  to  physical  processes  that  occur  in  nature:  a  strong  shock  wave  and  a  weak  shock 
wave  although  the  latter  is  preferred  by  nature  and  is  the  one  we  study  here.  The  third  solution  is  non¬ 
physical.  Frequently,  one  uses  a  numerical  root  finding  procedure  such  as  the  Secant  method  to  iteratively 
solve  Eq.(5.3)  for  the  value  of  (3  corresponding  to  the  weak  solution  but  here  we  used  the  analytic  solution 
(obtained  via  Mathematica  v4  [46]).  This  was  especially  useful  for  computing  exact  analytic  sensitivity 
derivatives  for  this  problem.  A  plot  of  the  (3  —  9  —  M  relationship  is  shown  on  the  left  of  Figure  5.2  for 
selected  Mach  numbers.  As  can  be  seen,  for  each  Mach  number  there  is  a  maximum  deflection  angle  beyond 
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Fig.  5.2.  (left)  (3  —  9  —  M  relationship  for  various  Mach  numbers,  (right)  pressure  rise  from  the  deterministic  solution. 


which  no  solutions  exist.  At  Mach  2,  the  maximum  deflection  is  approximately  23°,  at  Mach  5  it  is  roughly 
41°. 

Once  /3  is  found,  one  solves  Eq.(5.2)  for  the  normal  component  of  the  Mach  number  and  replaces  Mi  by 
Mi,^  in  Eq.(5.1)  to  obtain  the  deterministic  values  of  the  pressure  rise.  Other  than  the  assumptions  already 
mentioned  (i.e.,  inviscid,  perfect  gas)  which  arise  in  the  governing  equations,  this  procedure  is  exact.  The 
solution  for  the  pressure  rise  across  an  oblique  shock  is  shown  on  the  right  of  Eigure  5.2. 

A  far  more  general  albeit  expensive  approach  to  solving  this  problem  is  to  use  modern  CED  techniques. 
To  that  end,  we  used  the  two-dimensional  CED  code,  ELOW.f  (formerly  ANSERS)  supplied  by  A.  Taylor 
This  code  solves  the  Navier-Stokes  equations  and  subsets  including  the  Euler  equations  on  finite  volume 
meshes  using  upwind  discretization  methods  (see  [44] ) .  A  sample  result  of  one  calculation  is  shown  in  Eigure 
5.3  for  which  pressure  contours  over  the  wedge  corresponding  to  the  conditions  Moo  =  3  and  0  =  5°  are 
displayed.  To  compare  the  CED  calculations  with  oblique  shock  relations,  we  took  the  pressure  behind  the 
shock  to  be  the  cell-averaged  value  from  the  last  cell  on  the  surface  (i.e.,  right  before  the  exit  plane).  On  31 
X  46  grids,  the  results  of  the  CED  calculations  at  all  Mach  numbers  considered  are  visually  indistinguishable 
on  the  pressure  rise  plot  from  the  exact  oblique  shock  results. 

5.2.  Stochastic  Problem.  Two  variables,  the  Mach  number,  M,  and  wedge  angle,  9,  were  considered 
to  be  Gaussian  random  variables,  both  with  a  coefficient  of  variation  of  10%.  The  main  focus  was  on  the 
geometric  uncertainty  associated  with  the  wedge  angle.  The  uncertainty  analysis  methods  implemented 
were  Monte  Carlo  and  the  first-order  moment  method  for  which  the  first-order  sensitivity  derivatives  were 
computed  analytically.  Eor  example,  to  compute  the  sensitivity  derivative  ^  we  used  the  chain  rule, 

dp  dp  dMi^^  dj3 

^  M  ^  aMi„  dfd 

The  first  two  terms  are  trivial  but  the  third  is  essentially  impossible  by  hand.  Mathematica  version  4  [46]was 
used  to  evaluate  analytically,  although  other  approaches,  including  numerical  derivatives  could  be  used. 
The  results  of  the  sensitivity  derivative  evaluation  are  shown  in  Eigure  5.4.  Note  that  the  derivative  is  fairly 
linear  until  9  approaches  the  maximum  flow  deflection  angle  at  which  point  the  derivative  varies  very  rapidly. 
This  is  just  another  example  of  the  highly  non-linear  nature  of  fluid  mechanics. 
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Fig.  5.3.  N 011- dimensional  pressure  contours  for  Mach  S  flow  over  a  5°  wedge.  Second-order  accurate  upwind  solution 
from  FLOW.f  on  a  SI  x  46  mesh. 


Fig.  5.4.  Non-dimensional  sensitivity  derivative  of  pressure  with  respect  to  deflection  angle.  Note  the  highly  non-linear 
nature  in  the  vicinity  of  the  maximum  defection  angle  for  each  Mach  number.. 


Monte  Carlo  simulations  using  both  oblique  shock  relations  and  CFD  were  run  for  mean  Mach  numbers 
[M  =  2,3,  4,  5}  and  mean  flow  deflection  angles  {9  =  5°,  10°,  15°,  20°}.  For  each  case,  1000  simulations 
were  performed.  Additionally,  correlations  between  the  two  random  variables  were  investigated.  For  each 
{M.  9}  pair,  the  correlation  coefficient,  p,  was  varied  from  {—1, 1}  resulting  in  9  x  16  x  1000  =  144,000  2D 
CFD  simulations  and  an  equal  number  of  oblique  shock  calculations.  The  CFD  calculations  were  performed 
at  ICASE  on  8-processors  of  CORAL.  It  took  approximately  two  overnight  runs  to  complete  the  144, 000 
simulations.  The  oblique  shock  relations  were  run  in  a  matter  of  minutes  on  a  200  MHz  pentium  PC.  In 
terms  of  the  statistics  of  the  distributions,  both  approaches  yielded  consistently  similar  values.  The  expected 
value  of  the  pressure  rise  with  95%  mean  confidence  intervals  and  the  coefficient  of  variation  are  shown  in 
Figure  5.5. 

Note  that  the  data  point  corresponding  to  M  =  2,  9  =  20°  is  missing.  For  a  Gaussian  distributed 
random  variable  with  a  mean  of  20°  and  standard  deviation  of  2°  (i.e.,  10%  CoV),  the  95%  confidence  region 
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Fig.  5.5.  (left)  Mean  pressure  rise  for  an  input  CoV  [9]=0.1.  The  error  bars  represent  95%  mean  confidence  intervals 
from  1000  Monte  Carlo  simulations,  (right)  output  pressure  coefficient  of  variation.  The  “bad’'  point  at  Mach  2,  9  =  20^  is 
due  to  sampling  values  of  9  greater  than  9Max- 


Table  5.1 

Number  of  samples  from  N[20,2]  that  are  greater  than  9  (out  of  1000). 


9° 

Number  of  Samples  >  9 

Exact  Values 

Random  Draw 

20 

500 

492 

21 

309 

292 

22 

159 

143 

23 

67 

73 

24 

23 

24 

25 

6 

8 

26 

1 

2 

is  roughly  [16°,  24°].  Clearly,  when  randomly  drawing  1000  samples  from  such  a  distribution,  numerous 
samples  will  be  drawn  with  values  greater  than  23°  (see  Table  5.1)  which  corresponds  to  the  maximum  flow 
deflection  angle  for  this  Mach  number  hence  the  M  =  2,  9  =  20°  data  point  was  rejected.  To  further  illustrate 
this,  Figure  5.6  shows  the  input  histogram  with  mean  angle  of  20°  where  numerous  samples  greater  than  23° 
can  be  seen.  The  center  and  right  histograms  are  the  pressure  output  for  input  conditions  M  =  2,  6  =  15° 
and  M  =  5.  9  =  20°,  respectively  for  which  all  of  the  sampled  angles  are  well  below  9 Max  for  each  Mach 
number.  In  this  rather  academic  example,  we  chose  large  variations  in  the  input  variables  purposely  to  make 
the  graphics  easier  to  see,  but  for  one  case,  this  resulted  in  non-physical  sampling.  Such  is  likely  to  occur 
with  many  random  variables  that  may  arise  in  fluid  mechanics  applications,  for  example,  thermodynamic 
and  transport  properties  are  bounded  from  below  by  zero  whereas  a  Gaussian  distribution  has  support 
on  [— cx),cx)].  One  approach  to  overcoming  this  problem  is  to  either  truncate  the  distribution  or  to  use 
distributions  with  finite  support  (such  as  the  Beta  distribution)  provided  that  the  use  of  such  a  distribution 
can  be  justified. 

For  demonstration  purposes,  correlations  between  the  upstream  Mach  number  and  9  were  considered. 
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Fig.  5.6.  Histograms  of:  (left)  input  wedge  angles  from  N[20,2],  (center)  output  pressure  at  M  =  2,  9  =  15^  and  (right) 
output  pressure  at  M  =  5,  9  =  20*^. 


Fig.  5.7.  Output  pressure  coefficient  of  variation  due  correlation  of  free -stream  Mach  number  and  wedge  angle,  9. 

The  result  of  correlating  these  variables  is  shown  in  Figure  5.7.  As  expected,  when  M  and  0  are  positively 
correlated,  the  combined  model  uncertainty  increases  since  increasing  M  or  9  independently  increases  the 
uncertainty.  Likewise,  negative  correlation  reduces  the  variance  since  one  tends  to  offset  the  other.  Note  that 
a  rather  substantial  variation  with  correlation  coefficient  occurs  particularly  for  the  stronger  cases  (larger 
mean  M  and  mean  6).  Clearly,  when  considering  multiple  random  input  variables,  it  will  be  important 
to  have  knowledge  of  the  correlation  among  the  variables  in  order  to  get  reliable  output  estimates  of  the 
variance. 
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Expansion  Fan 


Fig.  6.1.  Sketch  of  supersonic  flow  over  an  expansion  corner. 

6.  Prandtl-Meyer  Expansion  Waves.  Expansion  waves  form  another  fundamental  building  block 
for  compressible  flow  theory.  In  this  case,  we  consider  the  flow  of  a  calorically  perfect  gas  across  a  centered 
expansion  fan  emanating  from  a  sharp  convex  corner  as  shown  in  the  sketch  of  Figure  6.1.  We  take  the 
output  quantity  of  interest  to  be  the  pressure  drop  through  the  expansion  fan  which  is  a  function  of  the 
upstream  Mach  number  and  flow  turning  angle. 

6.1.  Deterministic  Problem.  The  solution  to  this  supersonic  flow  problem  was  first  presented  by 
Prandtl  in  1907  and  subsequently  by  Meyer  in  1908  [5].  The  basic  problem  is  to  determine  properties  behind 
the  rearward  Mach  line  given  flow  properties  ahead  of  the  expansion  corner.  The  governing  differential 
equation  for  this  flow  is 

(6.1)  de  =  Vm2  - 1^ 

where  V  is  the  flow  velocity  and  d6  is  an  infinitesimally  small  expansion  angle.  Integrating  the  differential 
equation  over  the  entire  expansion  angle  yields, 

(6.2)  \M\  =  v{M2)  - 

where  the  notation  A0  implies  the  total  turning  angle  that  the  flow  experiences  and  v{M)  is  the  Prandtl- 
Meyer  function,  which  for  a  calorically  perfect  gas,  is 

(6.3)  viM)  =  \l^  ^  tan~^  ,1^  (M^  —  1)  —  tan~^  \/ —  1. 

V7— 1  V7+I 

Given  the  upstream  Mach  number.  Mi,  and  the  flow  expansion  angle,  |A0|,  the  procedure  for  solving  this 
problem  is: 

1.  compute  v{Mi)  from  Eq.(6.3) 

2.  compute  v{M2)  from  Eq.(6.2) 

3.  compute  M2  by  solving  Eq.(6.3)  as  a  root-finding  problem  or  find  M2  from  tabulated  values 

4.  Use  isentropic  relations  to  compute  flow  properties  in  region  2 
For  example,  the  pressure  decrease,  P2/P1,  can  be  obtained  from 

1+  2^M|1 

1  +  ^M2 


(6.4) 


P2_ 

Pi 
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Fig.  6.2.  N 011- dimensional  pressure  contours  for  Mach  S  flow  over  a  5^  expansion  corner.  Second-order  accurate  upwind 
solution  from  FLOW.f  on  a  SI  x  46  mesh. 


The  second  approach  that  we  used  to  solve  this  problem  was  the  flow  analysis  code  FLOW.f  described 
earlier.  Simple  h-meshes  consisting  of  31  x  46  grid  points  were  used  for  the  simulations.  Pressure  contours 
from  a  Mach  3  simulation  over  a  5°  expansion  corner  are  shown  in  Figure  6.2.  For  comparisons  with  the 
Prandtl-Meyer  results,  we  again  took  the  discrete  value  of  the  pressure  behind  the  expansion  to  be  the 
cell-averaged  value  at  the  last  cell  on  the  surface  ahead  of  the  exit  plane. 

6.2.  Stochastic  Problem.  Similar  to  the  oblique  shock  wave  example,  we  took  the  upstream  Mach 
number  and  the  flow  expansion  angle  to  be  Gaussian  random  variables  with  a  coefficient  of  variation  of  10%. 
Monte  Carlo  simulations  (1000  per  case)  and  the  first-order  moment  method  were  implemented  using  both 
analysis  methods  discussed.  The  expected  value  of  the  pressure  decrease  shown  in  Figure  6.3  closely  follows 
the  deterministic  solution.  On  the  original  21  x  31  grids  used  by  the  CFD  code,  we  found  discrepancies 
at  the  higher  Mach  numbers  and  turning  angles.  Consequently,  we  refined  the  grids  to  31  x  46  for  which 
solutions  at  all  [M,  9}  pairs  compared  closely  to  the  exact  solution. 

Sensitivity  derivatives  were  computed  using  finite-differencing  of  the  exact  inviscid  Prandtl-Meyer  solu¬ 
tion.  Results  for  the  pressure  sensitivity  are  shown  in  the  right  of  Figure  6.3.  For  moderate  turns  of  less 
than  approximately  10°,  pressure  sensitivity  increases  with  increasing  Mach  number,  but  at  higher  turning 
angles,  this  trend  reverses.  Since  the  first-order  estimate  of  the  standard  deviation  is  directly  proportional 
to  the  pressure  sensitivity  this  behavior  is  reflected  in  the  standard  deviation  of  the  pressure  distribution 
(see  Figure  6.4).  However,  the  relative  uncertainty  in  terms  of  the  coefficient  of  variation  monotonically 
increases  with  increasing  Mach  numbers  and  turning  angles.  Note  that  quite  substantial  variation  in  the 
output  occurs  at  the  higher  (M,  9}  pairs. 

We  implemented  the  bootstrap  estimate  of  the  standard  error  to  compute  confidence  intervals  for  any 
statistical  quantity.  The  method  has  the  advantage  that  it  is  easy  and  efficient  to  implement,  general  in  the 
sense  that  it  is  not  restricted  to  a  specific  distribution,  e.g.  a  Gaussian,  and  it  can  be  completely  automated 
for  any  estimator.  The  interested  reader  is  referred  to  [12]  for  details  of  the  method.  In  practice,  one 
takes  25  —  200  bootstrap  samples  to  obtain  a  standard  error  estimate.  Table  6.1  shows  bootstrap  standard 
error  estimates  of  the  mean,  standard  deviation,  and  coefficient  of  variation  for  a  range  of  bootstrap  sample 
numbers  for  the  case  of  Mach  5  flow  over  an  expansion  corner  with  a  mean  angle  of  20°  and  a  coefficient 
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Fig.  6.3.  (left)  Expected  value  of  the  pressure  decrease  through  an  expansion  from  1000  Monte  Carlo  simulations,  (right) 
Sensitivity  derivatives. 


Symbols  -  FLOW.f 


Fig.  6.4.  Pressure  output:  (left)  Standard  deviation  and  (right)  coefficient  of  variation  with  95%  mean  confidence  intervals 
from  1000  Monte  Carlo  simulations. 


of  variation  of  10%.  The  95%  confidence  intervals  shown  in  Figure  6.5  were  obtained  from  the  bootstrap 
estimate  of  the  standard  error  in  the  coefficient  of  variation  of  pressure. 

The  histogram  in  Figure  6.6  shows  that  the  random  variable  P2/P1  is  non-Gaussian  as  evidenced  by 
the  significant  tail  to  the  right.  Results  of  the  analysis  assuming  correlation  among  the  {M,6}  pairs  are 
presented  in  Figure  6.6.  The  trends  are  similar  to  the  oblique  shock  results  in  that  positively  correlated 
variables  result  in  the  greatest  uncertainty.  This  is  intuitive  since  increasing  the  expansion  angle  results  in 
greater  uncertainty  as  does  independently  increasing  the  Mach  number. 

7.  Supersonic  Airfoil.  Shock-expansion  theory  has  been  used  extensively  to  model  supersonic  flow 
over  thin  airfoils.  This  example  considers  steady,  two-dimensional  adiabatic  flow  of  an  inviscid,  perfect  gas 
over  the  airfoil  shape  described  by  its  upper  and  lower  surfaces,  yu{x)  and  yi{x),  respectively.  The  airfoil 
coordinates  are  given  by 

t  X  \ 

(7.1)  j/„(a;)  =  -cos[7r(- - yi{x)  = -yu{x) ,  a;G[0,c]. 
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Table  6.1 


Estimates  of  standard  error  for  various  bootstraps  samples  from  Mach  5  flow  over  an  expansion  corner  with  a  mean  angle 
of  20  degrees  and  a  coefficient  of  variation  of  10  percent. 


Number  of 

Bootstraps 

Standard  Error 

in  the  Mean 

Standard  Error 

in  (T 

Standard  Error 

in  CoV 

100 

0.000547105 

0.000751298 

0.0144571 

200 

0.000598759 

0.000673770 

0.0138239 

300 

0.000557259 

0.000678446 

0.0137223 

400 

0.000581007 

0.000674247 

0.0133587 

500 

0.000600424 

0.000691841 

0.0134268 

0.2  n 


0.15h 


>  0.1  h 

o 

o 


0.05  h 


Mach  5  Expansion  Flow 


Mean  CoV 
95%  Cl 


Q  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I 

5  6  7  8  9  10 

Flow  expansion  angle,  6 


Fig.  6.5.  Application  of  the  bootstrap  error  estimator  for  confidence  intervals. 


Fig.  6.6.  Left  -  Output  pressure  coefficient  of  variation  due  correlation  of  free-stream  Mach  number  and  wedge  angle,  9. 
Bight  -  Histogram  of  the  output  pressure  distribution  at  Mach  5  through  an  expansion  with  a  mean  angle  of  20^  and  a  10% 
coefficient  of  variation.  Note  that  the  distribution  is  skewed  to  the  right. 
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Fig.  7.1.  Surface  coordinates  (left)  and  pressure  distribution  corresponding  to  Moo  =  3  and  a  =  3°  (right). 


Fig.  7.2.  Grid  convergence  results  for  the  thin  supersonic  airfoil  at  Mach  3  and  an  angle  of  attach  of  3  degrees. 


The  shape  of  the  airfoil  is  shown  in  Figure  7.1.  Visually,  the  airfoil  is  displayed  five  times  thicker  than  it 
actually  is  due  to  the  aspect  ratio  of  the  plot.  For  all  calculations,  the  free-stream  Vlach  number  was  fixed  at 
three,  (Moo  =  3).  The  non-dimensional  pressure  distribution  obtained  on  a  uniform  65  grid  point  {h  =  0.016) 
mesh  at  an  angle  of  attack,  0  =  3°  and  a  thickness-to-chord  ratio,  i  =  0.05,  is  shown  on  the  right. 

Grid  convergence  studies  were  performed  prior  to  generating  stochastic  solutions.  Ten  grid  levels  were 
used  with  the  number  of  points,  n,  determined  from  n  =  2*+l,  i  =  {1,2,  ...10}.  The  mesh  size,  h  =  ^ 
varies  from  a  maximum  value  of  0.5  to  the  minimum  of  0.000976563.  Figure  7.2  shows  the  results  of  this 
grid  refinement  process  on  Ci,Cd,  and  Ci/Cd-  As  can  be  seen,  grid  convergence  occurs  fairly  quickly.  On  the 
65  grid  point  {h  «  0.016)  mesh  corresponding  to  the  sixth  level  of  refinement  (starting  from  a  very  coarse  3 
grid  point,  2  panel  airfoil),  the  solution  is  converged  for  all  practical  purposes. 

For  the  stochastic  calculations,  the  thickness-to-chord  ratio  was  assumed  to  behave  as  a  Gaussian  random 
variable  with  a  mean  value  of  ^  =  0.05  and  a  coefficient  of  variation  of  1%.  Solutions  on  a  65  grid  point  mesh 
were  generated  by  the  Monte  Garlo  method  and  the  first-  and  second-moment  methods.  The  sensitivity 
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x/c  Number  of  Simulations 


Fig.  7.3.  Monte  Carlo  convergence  history.  The  FOSM  and  SOSM  results  are  shown  for  reference. 


Fig.  7.4.  Expected  value  of  lift- to -drag  coefficients  with  uncertainty  due  to  a  1%  coefficient  of  variation  in  airfoil  thickness. 

derivatives  with  respect  to  ^  required  by  the  moment  methods  were  computed  using  second-order  accurate 
finite  differences.  The  upper  surface  pressure  sensitivity  derivatives  versus  x/c  for  two  different  angles  of 
attack  {a  =  0°  and  3°)  are  shown  in  Figure  7.3.  Three  function  evaluations  per  angle-of-attack  are  required 
by  this  procedure.  If  one  is  only  interested  in  a  single  output  quantity,  such  as  Ci/Cd,  there  is  no  need  to 
compute  the  pressure  sensitivities  since  the  sensitivity  derivative  of  the  output  quantity  can  be  determined 
directly  from  the  differencing  procedure.  Figure  7.3  shows  the  Monte  Carlo  convergence  of  the  expected 
value  of  CijCd  and  95%  mean  and  single  prediction  confidence  intervals. 

Finally,  the  performance  of  the  airfoil  across  a  range  of  angles  of  attack  was  studied.  The  angle-of- 
attack  was  varied  from  1°  to  5°  in  increments  of  1/4°.  Results  from  100  Monte  Carlo  simulations  for  each 
angle-of-attack  are  compared  with  results  from  the  first-order  moment  method  in  Figure  7.4  where  again 
good  correspondence  between  the  FOSM  confidence  interval  and  the  single  prediction  confidence  interval  is 
seen.  The  enlarged  view  on  the  right  clearly  shows  that  FOSM  slightly  underestimates  the  variance.  The 
second-order  moment  method  produces  results  virtually  identical  to  FOSM  and  were  omitted  for  clarity. 
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8.  Laminar  Boundary  Layer  Flow.  The  boundary  layer  equations  for  steady,  two-dimensional  in¬ 
compressible  flow  with  constant  properties  can  be  written  as 


(8.1) 


du  dv 
dx  ^  dy 


(8.2) 


du  du  dU  d‘^u 


subject  to  the  boundary  conditions 


(8.3) 


u{x,  0)  =  v{x,  0)  =  0  u{x,  cx))  =  U (x). 


It  is  well  known  that  these  partial  differential  equations  can  be  re-cast  in  terms  of  a  similarity  variable  as  an 
ordinary  differential  equation  for  which  solutions  exist  for  specific  freestream  velocity  distributions,  U{x). 
In  1908,  Blasius  found  a  solution  to  the  parallel  flow  over  a  flat  plate  with  U{x)  a  constant  (see  [45]).  Under 
these  conditions,  the  governing  equations  in  self-similar  form  reduce  to 

(8.4)  f"  +  ff"  =  0 
where  f  =  /(rj)  only  and 

(8.5)  u{x,y)  =  U{x)f'{ri) 

The  boundary  conditions  transform  to 

(8.6)  /(O)  =  /'(O)  =  0  /'(«))  =  1 


where  77(2;,  y)  =  y 


U 

2vx 


To  date,  an  analytic  solution  to  this  equation  is  not  known  although  series  solutions  with  a  finite  radius 
of  convergence  are  established.  Numerically,  the  problem  reduces  to  finding  the  value  of  /”(0)  such  that 
^  1  as  ?7  ^  CX).  An  asymptotic  analysis  shows  that  77  =  10  is  a  sufficiently  large  value.  Finding  the  correct 
value  of  is  a  simple  root  finding  problem  which,  in  our  calculations,  was  solved  by  the  Secant  method. 

The  correct  value  to  six  significant  digits  is  (see  [45]) 


(8.7)  /"(O)  =  0.469600  •  •  • 

8.1.  Deterministic  Problem.  A  relatively  fine  5000  grid  point  mesh  distributed  on  77  G  [0, 10]  was 
used  for  the  numerical  solution.  Similarity  functions  and  velocity  profiles  x  =  1  meter  for  three  different 
freestream  Reynolds  numbers  (Re  =  are  shown  in  Figure  8.1. 

8.2.  Stochastic  Problem.  The  kinematic  viscosity,  1^,  was  taken  as  a  normally  distributed  stochastic 
input  parameter  with  a  coefficient  of  variation  of  2%.  The  sensitivity  derivative  required  for  the  FOSM 
method  can  be  obtained  from 


(8.8) 


du  du  dr] 
dv  drj  dv 


The  second  derivative  (for  the  SOSAI  method)  can  be  found  by  further  application  of  the  chain  rule  to  be 


(8.9) 


0.0, 3/”, 
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Fig.  8.1.  Blasius  solution  for  the  flow  over  a  flat  plate.  Function  profiles  (left)  and  velocity  distributions  at  three  Reynolds 
tiumbers  (right). 


Fig.  8.2.  Non-dimensional  sensitivity  derivatives  (left)  and  distribution  of  the  velocity  standard  deviation  through  the 
boundary  layer. 


Plots  of  the  derivatives,  non-dimensionalized  by  U  and  the  mean  viscosity,  F,  are  shown  in  Figure  8.2. 
For  this  problem,  the  contribution  of  the  second  derivative  terms  to  the  estimates  of  the  mean  and  variance 
are  negligibly  small  and  hence  their  distributions  are  indistinguishable  as  seen  on  the  right  of  Figure  8.2 
where  the  Monte  Carlo  results  after  1000  simulations  are  also  shown  for  comparison. 

Mean  velocity  profiles  at  a;  =  1  meter  and  Re  =  10,  000  obtained  from  1000  Monte  Carlo  realizations 
are  shown  in  Figure  8.3.  An  enlarged  view  is  displayed  on  the  right  so  that  the  single  prediction  confidence 
intervals  can  be  readily  observed.  The  FOSM  and  SOSM  confidence  intervals  coincide  with  the  Monte 
Carlo  confidence  intervals  and  were  omitted  in  the  plot  for  clarity.  It  is  worth  noting  that  the  95%  mean 
confidence  intervals  are  smaller  by  a  factor  of  roughly  30  («  vlOOO)  and  should  be  used  for  giving  the  best 
estimate  of  the  mean.  However,  since  moment-based  methods  only  estimate  the  standard  deviation  of  the 
parent  population,  a  more  direct  comparison  to  Monte  Carlo  is  afforded  by  the  single  prediction  confidence 
intervals. 
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Fig.  8.3.  Velocity  profiles  through  the  boundary  layer  with  95%  single  prediction  confidence  intervals. 


Fig.  8.4.  Convergence  of  Monte  Carlo  simulations  at  y  —0.17  meters. 

The  convergence  of  the  standard  deviation  predicted  by  the  Monte  Carlo  method  versus  the  number 
of  simulations  is  shown  in  Figure  8.4.  The  FOSM  and  SOSM  estimates  (which  lay  on  top  of  one  another) 
are  also  shown.  Note  the  relatively  slow  convergence  of  the  Monte  Carlo  method  as  can  be  seen  by  the 
low  frequency  meandering  of  the  velocity  standard  deviation.  If  one  wants  very  accurate  estimates  of  the 
standard  deviation,  this  slow  convergence  and  hence  high  cost  will  be  unacceptable.  However,  if  one  only 
needs  a  moderately  accurate  estimate  of  the  mean  and  is  willing  to  accept  a  fairly  rough  estimate  of  the 
standard  deviation,  then  relatively  few  simulations  may  yield  acceptable  results. 

9.  Summary.  This  paper  reviews  deterministic  and  probabilistic  uncertainty  analysis  methods  applied 
to  fundamental  problems  in  fluid  mechanics.  Sources  of  uncertainty  and  a  discussion  of  selected  methods 
for  random  variable  models  are  presented.  Applications  presented  include  a  linear  convection  problem  with 
a  source  term,  a  non-linear  Burgers  equation,  supersonic  flow  over  wedges,  expansion  corners,  and  a  thin 
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supersonic  airfoil  as  well  as  incompressible  boundary  layer  flow. 

The  methods  discussed  and  implemented  are:  Interval  Analysis,  Propagation  of  error  using  sensitivity 
derivatives,  Monte  Carlo  simulation.  Moment  methods  and  Polynomial  Chaos.  Although  easy  to  imple¬ 
ment,  interval  analysis  often  results  in  maximal  error  bounds  that  are  quite  large.  The  basic  procedure  for 
implementing  Monte  Carlo  is  presented  next.  Although  computationally  intensive,  Monte  Carlo  solutions 
are  frequently  used  as  a  baseline  for  comparison  with  other  methods  since  they  are  known  to  converge  to 
the  exact  stochastic  solution  in  the  limit  of  infinite  sample  size.  First-  and  second-order  moment  methods, 
popular  because  of  the  relatively  low  cost  and  utility  in  a  design  environment  are  covered.  These  methods 
generally  yield  good  approximations  when  the  output  probability  density  function  is  a  Gaussian  distribution 
or  relatively  close  to  Gaussian.  Next,  Hermite  polynomial  chaos  is  described  for  solving  stochastic  problems 
involving  random  variables.  The  most  sophisticated  of  the  methods  reviewed,  polynomial  chaos  is  based  on  a 
spectral  representation  of  the  uncertainty  which  is  subsequently  decomposed  into  deterministic  and  random 
components.  Often,  highly  accurate  results  are  obtainable  from  this  approach  at  lower  cost  than  Monte 
Carlo  simulations. 

All  methods  were  implemented  for  a  non-linear  form  of  the  generalized  Burgers  equation  for  which  we 
obtained  an  exact  stochastic  solution.  To  mimic  the  behavior  of  CFD  codes,  we  used  second  order  spatial 
differencing  and  implemented  Newton’s  method  to  solve  the  non-linear  problem.  In  all  cases,  approximately 
ten  iterations  were  required  to  achieve  machine  precision  results.  Interval  analysis  error  bounds  were  unac¬ 
ceptably  large,  even  when  performed  as  a  single  function  evaluation  outside  the  iteration  loop.  Both  first- 
and  second-order  moment  methods  produced  reasonable  estimates  of  the  mean  and  variance  but  the  second 
order  estimates  were  substantially  better.  Polynomial  chaos  solutions  of  various  orders  were  generated.  The 
first-order  chaos  solutions  were  comparable  to  the  second  moment  solutions.  The  third  and  fourth-order 
solutions  were  very  accurate  and  matched  the  exact  PDF  of  the  solution  closely  at  all  points  in  the  domain. 
We  also  showed  that  the  treatment  of  boundary  conditions  and  the  quality  of  the  grid  has  an  impact  on  the 
error  convergence  as  a  function  of  the  order  of  the  chaos. 

Oblique  shocks  and  expansions  were  considered  in  which  the  flow  turning  angle  and  the  input  Mach 
number  were  considered  to  be  random  variables.  We  used  the  2  —  D  CFD  code  FLOW.f  as  well  as  shock 
expansion  theory  to  simulate  these  flows  by  the  Monte  Carlo  and  first-order  moment  methods.  Parametric 
results  showing  the  impact  of  correlation  among  the  variables  demonstrate  the  need  to  have  knowledge  about 
the  relationship  among  input  random  variables.  Supersonic  flow  over  a  thin  cosine  shaped  airfoil  was  studied 
at  a  variety  of  angles-of-attack.  The  thickness  of  the  airfoil  is  treated  as  a  random  variable  and  the  impact 
of  this  uncertainty  on  the  lift-to-drag  curve  is  examined. 

Finally,  we  studied  incompressible,  steady  flow  over  a  flat  plate,  the  celebrated  Blasius  flow,  in  which  the 
kinematic  viscosity  was  uncertain.  The  equations  were  solved  in  self-similar  variables  for  which  the  random 
variable  gets  folded  into  the  similarity  variable.  Both  moment  methods  and  Monte  Carlo  simulations  were 
performed.  For  this  example,  little  difference  was  observed  between  the  first-  and  second-moment  methods. 
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